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Abstract 

This  investigation  purports  to  develop  a  new  model  for  multiple  autonomous 
aircraft  mission  routing.  Previous  research  both  related  and  unrelated  to  this  en¬ 
deavor  have  used  classic  combinatoric  problems  as  models  for  Unmanned  Aerial 
Vehicle  (UAV)  routing  and  mission  planning.  This  document  presents  the  concept 
of  the  Swarm  Routing  Problem  (SRP)  as  a  new  combinatorics  problem  for  use  in 
modeling  UAV  swarm  routing,  developed  as  a  variant  of  the  Vehicle  Routing  Prob¬ 
lem  with  Time  Windows  (VRPTW).  The  SRP  removes  the  single  vehicle  per  target 
restraint  and  changes  the  customer  satisfaction  requirement  to  one  of  vehicle  on  lo¬ 
cation  volume.  The  impact  of  these  alterations  changes  the  vehicle  definitions  within 
the  problem  model  from  discrete  units  to  cooperative  members  within  a  swarm.  This 
represents  a  more  realistic  model  for  multi-agent  routing  as  a  real  world  mission  plan 
would  require  the  use  of  all  airborne  assets  across  multiple  targets,  without  con¬ 
straining  a  single  vehicle  to  a  single  target.  Solutions  to  the  SRP  problem  model 
result  in  route  assignments  per  vehicle  that  successfully  track  to  all  targets,  on  time, 
within  distance  constraints.  A  complexity  analysis  and  multi- objective  formulation  of 
the  VRPTW  indicates  the  necessity  of  a  stochastic  solution  approach  leading  to  the 
development  of  a  multi-objective  evolutionary  algorithm.  This  algorithm  design  is 
implemented  using  C++  and  an  evolutionary  algorithm  library  called  Open  Beagle. 
Benchmark  problems  applied  to  the  VRPTW  show  the  usefulness  of  this  solution  ap¬ 
proach.  A  full  problem  definition  of  the  SRP  as  well  as  a  multi-objective  formulation 
parallels  that  of  the  VRPTW  method.  Benchmark  problems  for  the  VRPTW  are 
modified  in  order  to  create  SRP  benchmarks.  These  solutions  show  the  SRP  solution 
is  comparable  or  better  than  the  same  VRPTW  solutions,  while  also  representing  a 
more  realistic  UAV  swarm  routing  solution. 
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Multi-Objective  UAV  Mission  Planning 


Using  Evolutionary  Computation 

I.  Introduction 

This  research  investigation  studies  the  routing  of  multiple  UAVs  and  UAV  swarms 
to  a  set  of  locations  while  meeting  constraints  of  time  on  target,  total  mission 
time,  enemy  radar  avoidance,  and  total  path  cost  optimization.  United  States  Air 
Force  (USAF)  research  has  been  focused,  and  increasingly  continues  to  focus,  on  the 
development  of  Autonomous  Unmanned  Aerial  Vehicles  (AUAV).  The  development 
of  UAV  technology  is  already  advanced,  with  uses  being  found  both  for  reconnaissance 
and  tactical  strikes  [52], 

While  hardware  for  these  UAV  systems  continues  to  evolve,  cutting  edge  re¬ 
search  focuses  on  the  development  of  taking  unmanned  systems  and  using  autonomous 
control  schemes  with  only  a  high  level  strategic  decision  maker  (human)  in  the  loop. 
Developing  a  single  autonomous  UAV  is  not  the  objective,  rather  the  objective  is  to 
develop  a  massive  array  of  AUAV,  capable  of  working  together  toward  a  common 
goal.  The  term  for  this  array  is  a  swarm  or  flock  for  which  there  are  many  differ¬ 
ent  design  approaches  as  the  problem  itself  exists  in  many  scientific  and  engineering 
domains.  Currently,  this  forefront  technology  exists  only  in  simulation  form  as  hard¬ 
ware  and  sensors  capable  of  implementing  an  actual  swarm  simply  do  not  exist.  This 
research  focuses  on  the  development  of  off-line  UAV  routing  and  mission  planning, 
combined  with  a  simulation  and  visualization  system  the  purpose  of  which  is  to  better 
understand  the  computational  complexity  of  multi  AUAV  routing. 

1.1  Problem  Statement 

The  problem  of  mission  planning  consists  of  assigning  multiple  vehicles  sets  of 
targets  to  visit.  These  targets  exist  in  a  field  of  uneven  terrain  where  different  enemy 
radar  lines  of  sight  exist.  There  exist  two  problem  aspects  to  deal  with,  the  first  is 
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the  development  of  flight  paths  between  targets.  The  path  must  be  optimized  for 
cost  and  risk.  Cost  is  how  much  energy  or  time  it  takes  to  traverse  the  path  and  risk 
is  a  measure  of  how  dangerous  the  flight  area  is.  The  second  is  the  development  of 
path  order.  Once  it  is  determined  how  to  best  fly  between  targets  the  order  of  these 
flight  paths  must  be  determined.  Generating  the  cost  of  the  path  is  a  separate  and 
immaterial  problem  related  to  the  development  of  path  order.  In  fact,  the  development 
of  single  path  optimization  is  already  being  studied  in  fields  such  as  robotics,  land, 
and  air  based  agents.  Once  these  path  costs  are  known,  or  estimated,  however  what 
development  process  should  be  used  to  structure  a  valid  route  plan  from  them?  This 
process  of  route  development  is  the  subject  matter  of  this  investigation. 

In  order  to  model  this  routing  situation,  a  combinatorics  problem  known  as  the 
Vehicle  Routing  Problem  with  Time  Windows  (VRPTW)  is  used.  The  VRPTW  en¬ 
compass  a  situational  problem  composed  of  a  number  of  vehicles,  known  targets  with 
time  visitation  constraints,  and  constraint  on  the  visitation  capacity  of  each  vehicle. 
This  model  most  efficiently  possesses  all  the  aspects  of  the  problem  under  considera¬ 
tion  and  is  well  documented  and  understood.  The  VRPTW  is  limited  in  its  ability  to 
model  realistic  UAV  routing,  necessitating  in  its  extension  into  a  new  problem  model. 
This  VRPTW  variation  is  called  the  Swarm  Routing  Problem  (SRP),  presented  as  a 
more  efficient  form  to  model  the  routing  of  multiple  UAVs  to  multiple  targets  con¬ 
currently.  This  problem  structure  is  illustrated  in  Figure  1.1.  Note  that  the  path 
planner  supplies  cost  values  to  the  router  portion  of  the  structure.  What  these  values 
are  is  irrelevant  to  the  router  itself.  The  purpose  of  the  router  is  to  efficiently  search 
the  innumerable  combinations  of  costs  supplied  to  in  order  to  determine  the  most 
efficient  (or  cost  effective)  path  possible.  The  remainder  of  this  section  introduces 
these  problem  concepts  in  greater  detail. 

1.1.1  The  Vehicle  Routing  Problem  with  Time  Windows.  The  Vehicle  Rout¬ 
ing  Problem  (VRP)  is  a  classic  combinatorics  problem  in  the  field  of  computer  science 
and  the  VRPTW  is  a  variant  of  it.  In  this  context  it  models  the  routing  of  UAVs 
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Figure  1.1:  Problem  Structure 


to  a  set  of  targets  within  a  time  window  for  each  target  and  an  overall  time  limit 
for  returning  to  the  start  location.  The  objective  of  the  problem  is  to  develop  a  set 
of  routes  for  a  fleet  of  vehicles  where  each  vehicle  services  a  target  within  its  time 
window  and  returns  to  the  start  point  after  all  the  targets  in  its  route  have  been 
serviced.  The  problem  is  constrained  by  how  many  targets  a  vehicle  can  visit,  each 
target’s  time  window,  and  a  return  time  limit  to  the  start  location  [51]. 

1.1.2  The  Swarm  Routing  Problem.  The  SRP  is  a  variation  on  the  VRPTW 
defined  for  the  first  time  in  this  work.  While  the  classic  VRPTW  restricts  each  vehicle 
to  a  single  route  and  target,  the  SRP  assumes  the  vehicles  can  join  together  and  split 
apart  at  different  locations,  in  essence  treating  each  vehicle  as  a  member  of  a  sub 
swarm  as  opposed  to  a  discrete  vehicle.  The  effect  of  using  the  SRP  to  model  UAV 
swarm  mission  routing  is  that  a  more  efficient  and  realistic  deployment  plan  is  created. 
Basing  the  model  on  the  VRPTW  allows  for  the  use  of  previously  made  insights  into 
the  VRP  and  VRPTW  to  create  effective  solutions  for  the  SRP.  Analysis  of  the 
VRPTW  also  allows  for  the  development  of  comparable  results  during  the  solution 
evaluation  phase,  the  result  of  which  is  a  metric  of  performance  to  compare  the  SRP 
results  to.  The  SRP  problem  complexity  is  NP-Hard  by  merit  of  its  derivation  from 
the  VRPTW.  Chapter  111  contains  a  more  exact  complexity  analysis. 
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1.2  Routing  Model  Assumptions 

The  routing  model  (for  both  VRPTW  and  SRP)  is  viewed  as  a  static  combi¬ 
natorics  problem.  While  the  VRPTW  model  does  take  into  account  time  constraints 
for  targets  that  is  also  all  it  is  capable  of  doing.  Modeling  the  routing  problem  as  a 
VRPTW  imposes  the  limitation  that  the  solution  only  be  for  the  VRPTW.  As  such, 
no  dynamic  aspects  can  be  introduced  such  as  in-mission  rerouting  or  defensive  mea¬ 
sures.  The  problem  does  not  take  possible  losses  of  aircraft  during  the  mission  into 
account.  This  limits  the  usefulness  of  the  model  as  a  real  world  problem  solver,  but 
is  a  closer  approximation  than  a  more  simplistic  point  to  point  router  without  time 
constraints. 

The  VRPTW  and  SRP  are  also  generalized  combinatorics  problems  and  no 
attempt  is  made  integrate  real  world  problem  constraints  such  as  vehicle  movement 
constraints,  real  world  distance  constraints,  or  communication  issues.  The  solution  to 
these  general  models  represents  a  top  level  routing  strategy  where  the  actual  values 
of  individual  route  sections,  or  how  they  are  determined,  is  immaterial.  In  reference 
to  the  routing  problem,  there  are  only  the  individual  path  cost  values,  what  they 
actually  are  or  how  they  are  developed  exists  only  within  the  users  perspective. 

1.3  Research  Concepts 

The  research  concepts  developed  as  part  of  this  investigation  involved  the  fol¬ 
lowing  solution  approaches.  They  are  introduced  to  give  the  reader  a  better  under¬ 
standing  of  the  direction  taken  in  this  endeavor. 

1.3.1  Multi- objective  Evolutionary  Algorithms.  Due  to  the  high  complexity 
of  the  problem  domain  in  this  investigation,  an  evolutionary  algorithm  approach  is 
used  as  previous  work  that  has  shown  a  great  deal  of  success  through  its  utilization 
[1,33,40,44,60].  Evolutionary  algorithms  encompass  a  number  of  different  stochastic 
search  techniques  where  random  solutions  are  “bred”  into  more  effective  solutions 
through  a  process  of  selection  and  mutation.  This  work  examines  how  multi-objective 
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approaches  obtain  the  most  effective  results  for  highly  constrained  solution  spaces, 
such  as  those  dealt  with  in  this  research.  The  concept  of  multi-objective  evolutionary 
search,  non-dominated  (Pareto)  front  solutions,  as  well  as  alternative  methods  of 
multi-objective  optimization  is  further  examined  in  Chapter  II. 

1.4  Research  Goals,  Objectives,  and  Assumptions 

The  goal  of  this  research  centers  on  the  creation  of  new  routing  software  for 
the  UAV  swarm  routing  problem  developed  in  previous  research  [11,44],  The  solu¬ 
tion  is  designed  in  the  context  the  it  integrate  into  previous  work  in  path  planning 
development  and  UAV  simulation. 

1.4-1  Objective  1:  Develop  SRP  Multi- objective  Problem  Model.  The  first 
objective  is  to  develop  a  problem  model  for  multiple  UAV  routing.  The  VRPTW 
serves  as  a  basis  for  the  development  of  this  model.  The  purpose  in  this  is  to  create  a 
problem  that  serves  as  a  more  realistic  model  for  developing  time  constrained  routes 
to  multiple  targets.  The  problem  is  formulated  in  terms  of  the  multiple  objectives  of 
total  path  length,  vehicle  count,  vehicle  wait  time,  and  average  path  length. 

1-4-2  Objective  2:  Develop  and  validate  MOEA  solution  to  VRPTW.  The 
second  objective  is  to  develop  a  multi-objective  evolutionary  solution  to  the  VRPTW. 
This  is  accomplished  by  examining  previous  MOEA  approaches  to  the  VRP  and 
VRPTW  in  order  to  develop  effective  genetic  operators  for  a  new  implementation. 
Completion  of  this  objective  consists  algorithm  development  and  fully  defining  genetic 
operators  the  Genetic  Algorithm  (GA)  uses.  The  purpose  in  developing  a  solution 
to  the  VRPTW  is  to  obtain  a  reference  point  for  development  of  the  SRP.  The  SRP 
represents  an  untested  model,  therefore  the  most  intelligent  action  is  to  validate  the 
solution  procedure  on  a  pre-existing  problems. 

A  variety  of  different  benchmark  problems  are  solved  using  the  implemented 
algorithm  and  the  results  are  compared  to  those  found  in  the  literature.  The  com- 
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plction  of  this  objective  is  signified  by  the  MOEA  approach  achieving  as  good  or 
better  results  than  those  found  in  the  literature  and  a  documentable  distribution  of 
solutions  across  the  objective  landscape.  The  objectives  of  comparison  are  total  path 
length,  vehicles  used,  and  total  wait  time  for  all  vehicles.  Each  of  these  objectives 
then  represent  metrics  of  performance  meant  to  be  minimized.  Experimental  analysis 
consists  of  a  variety  of  algorithmic  approaches  in  order  to  examine  the  utility  of  the 
each  method  and  statistical  analysis  is  used  to  show  the  significance  of  the  different 
methods  used. 

1.4-3  Objective  4-'  Develop  and  validate  MOEA  solution  to  SRP.  The  forth 
objective  is  to  develop  a  solution  to  the  SRP.  This  new  problem  is  a  variation  of  the 
VRPTW  and  is  solved  in  a  similar  fashion.  A  set  of  benchmark  problems  are  created 
by  altering  existing  VRPTW  problems.  These  benchmark  problems  are  then  solved 
using  a  solution  designed  parallel  to  the  one  developed  for  the  VRPTW. 

1.5  UAV  Model  Assumptions 

No  attempt  is  made  to  completely  model  the  aerodynamics  of  the  UAVs  used. 
They  are  assumed  to  be  point  masses  whose  control  scheme  is  some  undefined  method, 
though  the  scheme  would  most  likely  be  some  form  of  self-organized  control.  All 
routing  operations  take  place  off-line.  It  is  assumed  that  the  UAVs  have  some  form 
of  communication  method  though  no  attempt  is  made  in  this  work  to  model  the  type 
of  communication  architecture  the  swarm  may  use.  It  is  assumed  control  and  sensor 
mechanisms  exist  that  allow  the  UAVs  to  pilot  themselves  without  the  need  for  human 
management. 

References  throughout  are  made  to  UAVs.  The  term  is  meant  to  refer  to  au¬ 
tonomous  UAVs  (AUAV),  not  remote  controlled  vehicles.  Chapter  II  contains  back¬ 
ground  information  on  self  organized  control  theory  which  is  presented  as  a  possible 
control  method  for  AUAVs  traversing  through  a  routing  solution. 
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1.6  Sponsor 

This  research  is  sponsored  by  the  Air  Force  Research  Laboratory  (AFRL)  Sen¬ 
sors  Applications  and  Demonstrations  Division  (AFRL/SNZ),  specifically,  the  Virtual 
Combat  Laboratory  (AFRL/SNZW)  at  Wright  Patterson  Air  Force  Base,  Ohio.  The 
Virtual  Combat  Laboratory  conducts  advanced  development  initiative  as  well  as  field 
and  flight  test  demonstrations  and  evaluations.  A  series  of  UAV  simulation  models 
are  maintained  within  the  lab  along  with  a  suite  of  visualization  tools.  This  research 
continues  the  ongoing  relationship  between  AFIT/ENG  and  AFRL/SNZW  in  which 
both  parties  share  information  on,  and  enhance  the  capabilities  of,  UAV  swarm  sim¬ 
ulations.  The  point  of  contact  at  AFRL  is  Mr.  Mike  Foster  (AFRL/SNZW). 

1.7  Thesis  Construction 

This  thesis  document  has  been  constructed  with  the  goal  that  it  be  readable, 
efficient  and  thorough.  The  basic  concepts  of  research  are  listed  in  this  chapter  and 
expanded  upon  in  Chapter  II.  A  thorough  examination  of  the  research  problem  under 
consideration  is  presented  in  Chapter  III.  High  level  solution  composition  is  detailed 
in  Chapter  IV.  The  implementation  documentation  is  contained  in  Chapter  V  and 
Appendix  B.  Chapter  VI  details  the  experiments  used  to  validate  the  solution  tech¬ 
niques  proposed.  Chapter  VII  contains  results  and  analysis.  Chapter  VIII  concludes 
the  thesis  with  a  discussion  of  contributions  suggested  future  endeavors. 
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II.  Background  Research 


This  chapter  is  divided  into  several  sections  relevant  to  the  areas  of  research  pur¬ 
sued.  The  first  section  discusses  previous  work  that  this  research  extends.  Each 
subsequent  section  is  presented  as  an  examination  of  relevant  background  information. 
The  background  information  also  contains  problem  domain  formulations  of  the  vehicle 
routing  and  path  planning  problems.  The  second  section  is  an  analysis  of  the  Vehicle 
Routing  Problem  with  Time  Windows  (VRPTW)  and  how  its  solution  can  be  used 
to  effectively  solve  the  Unmanned  Aerial  Vehicle  (UAV)  routing  problem  for  multiple 
targets.  The  third  section  is  an  examination  of  the  path  following  constraint  that 
is  formulated  in  conjunction  with  the  VRPTW.  The  fourth  section  covers  the  topic 
of  Multi-objective  Evolutionary  Algorithm  (MOEA)s.  Pareto  ranking  techniques  are 
discussed  as  well  as  the  benefits  of  using  MOEA  approaches  over  weighted  and  single 
objective  techniques.  The  chapter  concludes  with  a  section  on  the  concept  of  self 
organized  UAV  swarms.  Previous  work  in  this  area  as  well  as  its  context  in  this  work 
is  examined. 

2.1  AUAV  Simulation  Research  History  at  AFIT 

Research  into  the  problem  of  effectively  routing  and  controlling  a  swarm  of 
UAVs  has  been  under  active  development  since  2000  (at  AFIT).  The  problem  is  by 
no  means  simple  and  crosses  into  many  scientific  fields.  Currently,  the  technology 
exists  mainly  in  the  simulation  arena,  though  much  work  has  been  done  with  single 
UAV  control  [1].  The  following  discussion  introduces  previous  thesis  research  efforts 
that  have  dealt  with  the  topic  of  UAV  routing  and  simulation.  This  information  is 
summarized  in  Table  2.1. 

Research  originally  began  with  UAV  routing  by  modeling  the  problem  as  a  clas¬ 
sic  combinatoric  problem  called  the  Traveling  Salesman  Problem  (TSP).  In  Sezer  [43], 
an  advanced  search  tree  method  was  developed  for  solving  a  constrained  mission  plan¬ 
ning  problem.  Research  then  moved  into  stochastic  methods  when  Secrest  [42]  devel¬ 
oped  a  particle  swarm  optimization  method  for  solving  the  TSP.  At  this  point  it  had 


Table  2.1:  Previous  AFIT  theses  related  to  AUAV  routing. 


Title 

Author 

Year 

Subject 

Mission  Route  Planning  with 

Multiple  Aircraft  &  Targets  Using 

Parallel  A*  Algorithm 

Ergin  Sezer 

2000 

Using  a  tree  search  method  to  solve  routing 

problems 

Traveling  Salesman  Problem  for 

Surveillance  Mission  Using  Particle 

Swarm  Optimization 

Barry  Secrest 

2001 

Using  a  combinatorics  problem  as  a  model  for 

UAV  routing 

Multi-objective  Mission  Route 

Planning  Using  Particle  Swarm 

Optimization 

Ekursat  Yavuz 

2002 

Using  PSO  search  method  to  route  UAVs 

Distributed  Control  of  a  Swarm  of 

Autonomous  Unmanned  Aerial 

Vehicles 

James  Lotspeich 

2003 

Examining  the  control  problem  of  UAV 

swarms 

Swarming  Reconnaissance  Using 

Unmanned  Aerial  Vehicles  In  A 

Parallel  Discrete  Event  Simulation 

Joshua  Corner 

2004 

Developing  the  AFIT  swarm  simulator 

Evaluation  and  Optimization  of 

UAV  Swarm  Multi-Objectives 

Mark  Kleeman 

2004 

Developing  the  multi-objective  idea  of  UAV 

routing 

Evolution  of  Control  Programs  for 

a  Swarm  of  Autonomous 

Unmanned  Aerial  Vehicles 

Keven  Milam 

2004 

Using  GP  to  develop  the  control  structure 

A  Genetic  Algorithm  for 

Unmanned  Aerial  Vehicle  Routing 

Matt  Russell 

2005 

Applied  new  chromosome  structure  to  the 

VRP 

Evolving  Self  Organizing  Behavior 

for  Homogeneous  and 

Heterogeneous  Swarms  of  UAVs 

and  UCAVs 

Ian  C.  Price 

2006 

Optimizing  the  control  structure  of  UAVs 

AFIT  UAV  Swarm  Mission 

Planning  and  Simulation  System 

James  Slear 

2006 

Development  of  path  planning  software 

been  established  that  the  high  complexity  of  the  routing  problem  necessitated  the  use 
of  stochastic  methods  in  order  to  achieve  results  approaching  optimality  in  a  tractable 
time  frame.  In  Yavuz  [59],  particle  swarm  optimization  was  used  to  develop  optimal 
solutions  to  a  multi-objective  routing  problem.  Based  on  lessons  learned  in  that  thesis 
effort  (of  the  different  objectives  relevant  to  the  problem),  further  work  was  done  by 
Russell  [40]  where  the  problem  was  formulated  as  a  VRP  and  solved  using  a  GA  [7]. 
Slcar  [44]  used  this  problem  formulation  while  path  planning  functionality  was  added. 
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In  Kadrovach  [23],  Lotspeich  [28],  and  Milam  [31]  methods  for  Self  Organized 
(SO)  control  were  developed.  These  control  schemes  were  based  on  flock  and  swarm 
behaviors  first  defined  by  Reynolds  [38] .  Simulations  using  these  control  schemes  were 
at  first,  two  dimensional  but  illustrated  the  validity  of  the  idea  behind  flocking  and 
swarming. 

After  it  was  shown  that  self  organized  control  could  be  effectively  implemented, 
research  diverged  into  three  areas:  swarm  routing,  swarm  communication,  and  SO  al¬ 
gorithm  development.  Swarm  communication  was  approached  by  Kadrovach  [23]  and 
Kleeman  [24]  where  it  was  shown  that  adequate  communication  within  a  real  world 
swarm  is  an  extremely  complex  problem.  Disregarding  the  hardware  limitations  of 
current  wireless  communication  technology  and  the  number  of  external  sensors  that 
would  be  required,  simply  developing  the  algorithms  that  would  be  able  to  handle  the 
massive  amount  of  internal  swarm  communication  is  still  a  largely  unsolved  problem. 
Current  developments  in  self  organized  control  have  centered  on  creating  MOEAs 
which  optimize  the  SO  rule  settings.  Price  [36]  developed  an  evolutionary  algorithm 
and  simulator  that  was  able  to  optimize  attack  patterns  of  the  swarm  in  a  2D  simu¬ 
lation.  Different  attack  patterns  as  well  new  optimization  strategies  were  created  by 
Nowak  [32], 

Both  Russell  [40]  and  Slear  [44]  utilized  the  AFIT  UAV  simulator  which  was 
first  developed  by  Kadrovach  [23]  and  then  ported  into  a  Linux  environment  and 
parallelized  in  Corner  [11],  Corner  [11]  also  developed  a  Parallel  Discrete  Event 
Simulation  (PDES)  methodology  to  simulate  larger  and  more  complex  SO  swarms  in 
3D. 

Many  different  languages  and  methodologies  have  been  applied  to  the  develop¬ 
ment  of  routing,  control  and  simulation  solutions.  The  most  popular  coding  languages 
(in  related  thesis  work  and  in  published  research)  have  been  Matlab,  Java™,  and 
C++.  A  2D  flocking  simulation  in  Matlab  was  created  by  Watson  [54]  to  study  the 
dynamics  of  the  self  organizing  behavior.  A  similar  simulation  developed  in  Java  was 
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created  by  Price  [36]  which  treated  the  UAVs  as  point  masses  with  no  dynamic  control 
mechanisms.  Previous  routing  [40]  and  path  planning  programs  [44]  have  been  cre¬ 
ated  using  C  and  C++  and  both  utilized  an  evolutionary  computation  library  called 
GALib  [53]  (also  written  in  C++),  as  the  support  library  for  their  solution  techniques. 
The  simulator  software  is  based  on  a  C++  library  called  SPEEDES,  which  was  de¬ 
veloped  to  provide  communication  support  for  PDES.  The  selection  of  this  library 
was  based  on  an  investigation  completed  by  Corner  [11], 

The  problem  of  what  a  swarm  actually  does  when  it  reaches  a  target  (i.e.  mission 
type)  has  only  recently  been  dealt  with  due  to  the  routing  problem  itself  being  such 
a  complex  endeavor.  Some  recent  work  has  been  conducted  in  2D  environments  that 
simulate  different  UAV  missions  [36].  Currently,  UAV  missions  are  viewed  in  a  purely 
reconnaissance  oriented  role.  With  the  capabilities  of  a  real-world  swarm  there  are 
no  limits  to  what  its  capabilities  would  be.  A  swarm  could  be  sent  out  to  perform 
a  concentrated  attack  on  a  target,  patrol  for  and  attack  targets  of  opportunity,  or 
even  escort  manned  targets  across  hostile  areas  by  not  only  targeting  enemy  units 
but  patrolling  the  area  to  alter  the  course  of  the  manned  unit  with  real  time  data. 

The  approach  to  UAV  swarm  simulation  has  been  very  iterative  in  development, 
first  with  problem  development  and  understanding,  then  routing  and  control  solutions, 
and  more  recently  advanced  simulations.  This  research  also  serves  as  an  iterative  step 
to  improve  on  the  problem  model  used  to  develop  routing  solutions,  apply  true  multi¬ 
objective  techniques  to  the  solving  routing  and  path  planning  problem  and  enhance 
the  current  simulator  technology.  The  basic  problem  remains  the  same  now  as  it  has 
always  been: 

Develop  a  method  for  taking  a  set  of  targets  within  an  area  of  terrain  and  finding 
a  set  of  paths  for  each  UAV  or  UAV  swarm  such  that  all  sites  are  visited  by  an 
appropriate  number  of  UAVs  to  accomplish  the  desired  mission  while  minimizing  path 
cost  and  maximizing  UAV  safety. 
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2.2  Vehicle  Routing  Problem  with  Time  Windows 

The  VRP  is  a  well  established  combinatorics  problem  with  many  variations  and 
solutions  [51],  one  of  these  variations  being  the  Vehicle  Routing  Problem  with  Time 
Windows  (VRPTW)  [51].  The  problem  is  most  easily  explained  with  a  real  world 
example.  There  exists  a  package  delivery  service  warehouse  that  has  a  number  of 
trucks.  Each  truck  can  carry  a  limited  load  of  packages  to  customers.  Customers  are 
spread  out  around  the  depot  and  each  can  only  be  visited  during  a  certain  time  of  the 
day.  Likewise,  the  warehouse  itself  is  only  open  for  a  certain  number  of  hours  until 
it  closes.  Determine  the  shortest  set  of  paths  for  each  vehicle  to  take,  such  that  all 
customers  are  visited  within  their  window  of  opportunity  and  all  vehicles  return  on 
time. 

From  this  example,  it  is  seen  that  the  VRPTW  consists  of  a  set  of  targets,  some 
number  of  vehicles,  and  a  depot.  The  depot  is  the  deployment  and  return  point  for 
all  the  vehicles.  Each  target  (and  the  depot)  has  a  Euclidian  location  (coordinate), 
some  associated  demand  (except  the  depot),  an  arrival  time  window,  and  a  path  to 
every  other  location.  The  objective  of  the  problem  is  to  develop  a  set  of  routes  for 
each  vehicle,  so  that  all  targets  are  visited  within  the  time  window,  the  associated 
demand  is  met,  and  all  vehicles  return  to  the  depot  on  time.  The  solved  problem 
is  visualized  in  Figure  2.1  where  the  darker  areas  correspond  to  the  time  window  in 
which  a  vehicle  can  visit  the  target.  Each  vehicle  in  the  problem  has  a  set  capacity 
that  it  can  not  exceed  while  visiting  customers.  Visiting  a  customer  subtracts  its 
demand  from  the  vehicles  capacity.  The  total  demand  on  the  vehicle  is  the  sum  of 
the  demands  of  all  the  customers  visited.  If  the  vehicles  being  used  do  not  have 
some  type  of  capacity  constraint  the  problem  then  decomposes  to  a  TSP  since  one 
vehicle  can  now  satisfy  all  customers.  The  remainder  of  this  section  is  a  more  detailed 
mathematical  description  of  the  single-objective  VRPTW  presented  in  order  to  better 
understand  the  constraints  and  optimization  goal  of  the  problem. 


12 


Figure  2.1:  The  Vehicle  Routing  Problem  with  Time  Windows. 

The  VRPTW,  as  formulated  by  Toth  [51]  based  on  the  ordinal  formulation 
by  Solomon  [45],  is  defined  by  a  fully  connected  graph  G  =  (V,  A)  where  V  = 
{v0,vi,  ...,vn}  and  v0  is  the  depot.  The  set  of  edges  is  defined  as  A  =  {( Vi,Vj )  G 
V,i  t -  j}  where  each  edge  has  associated  with  it  some  travel  cost  c(vi,Vj).  The  cost 
of  the  edge  is  the  distance  from  target  i  to  target  j  in  the  context  of  the  problem 
discussed  here.  Cost  is  not  necessarily  distance  however,  cost  is  only  a  value  that  is 
associated  with  traversing  the  route.  It  could  be  the  distance,  it  could  be  the  time, 
or  it  could  be  the  amount  of  fuel  required  to  travel  the  path.  A  time  window  exists 
for  all  customers  where  E  represents  the  earliest  start  time  and,  L ,  the  latest  arrival 
time.  The  latest  arrival  time  is  the  point  at  which  the  vehicle  can  arrive  and  still 
complete  the  service  time  defined  by  S.  If  the  vehicle  arrives  earlier  than  A,  it  incurs 
a  waiting  time,  W,  which  is  the  difference  between  the  arrival  time  and  E.  The  total 
time  a  vehicle  takes  to  complete  it’s  route  is  the  summation  of  all  route  path  travel 
costs,  waiting  times,  and  service  times  Wi  +  si)  and  is  refereed  to  as 

the  cost  of  the  path.  The  total  path  cost  must  not  exceed  the  latest  arrival  time  (i.e. 
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closing  time)  of  the  depot.  Cost  is,  in  essence,  the  total  time  required  to  traverse  the 
route. 

Each  customer  also  has  associated  with  it  some  demand,  D,  satisfied  by  the 
vehicle  servicing  it.  There  exists  a  number  of  vehicles,  K,  where  each  vehicle  has 
some  capacity  constraint,  C,  such  that  the  sum  of  the  demands,  D,  for  all  customers 
visited  by,  k  G  K,  does  not  exceed  C.  The  solution  to  the  problem  is  a  list  of  ordered 
targets  for  each  vehicle,  such  that  the  visitation  of  each  target  fulfills  all  customer 
needs  without  violating  any  time  or  demand  constraints.  The  objective  then  is  to 
determine  the  set  of  paths  for  the  vehicles  such  that  the  total  cost  is  minimized.  Note 
that  the  cost  of  a  route  is  not  the  total  time  the  route  takes  to  complete,  it  is  only 
the  sum  cost  of  the  edges  the  vehicle  traverses. 

The  following  is  a  mathematical  programming  formulation  of  the  single-objective 
VRPTW,  again  taken  from  [51].  Vehicles  are  defined  within  the  problem  by  their  in¬ 
clusion  in  a  flow  variable  which  is  a  binary  value  indicating  if  vehicle  k  exists  on 
the  path  that  connects  (i,j)  G  V  at  any  point  in  the  solution.  A  time  variable 
indicates  the  start  time  of  vehicle  k  at  location  i.  The  subscript  j  G  A ±(z)  indicates 
the  set  of  edges  from  i  to  j  where  j  is  not  equal  to  i,  the  plus  or  minus  indicates  either 
a  forward  or  backward  move  along  the  path. 

Aij  -  Edge  cost  between  i  and  j 

Vn  -  Network  vertices  for  n  customer  (Vo  is  the  depot) 

En  -  Earliest  arrival  time  of  customer  n 
Ln  -  Latest  arrival  time  of  customer  n 
Sn  -  Service  time  of  customer  n 
Dn  -  Demand  of  customer  n 
K  -  Set  of  Vehicles 
Ck  -  Capacity  of  vehicle  k 
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Xijk  ^  {0,  1} 


Vfc  g  AT,  ( i,j )  g  A, 


(2.1) 


a;ifc  >  0  ieN,keK ,  (2.2) 

Equations  (2.1)  and  (2.2)  define  the  flow  and  time  variables  used.  The  flow 
variable  x^k  is  a  binary  value  that  indicates  vehicle,  k,  travels  from  location,  i  to  j, 
if  equal  to  one,  and  zero  otherwise.  The  time  variable,  u;^,  specifies  the  start  time  at 
location,  i,  by  vehicle  k. 


E  E  x^k  =  1  Vi  G  N,  (2.3) 

k&K  je A+(i) 

y,  Xojfc  =  1  Vfc  G  K,  (2.4) 

jeA+(o) 

y  Zi,n+i,fc  =  1  Vfc  G  A',  (2.5) 

ieA~(n+l) 

y  Xjjfc-  y  x^k  —  0  VA;  G  AT,  Vi  G  N,  (2.6) 

*GA+(i)  ieA~(j) 


Equations  (2.3)-(2.6)  define  the  edge  constraints  of  the  graph  in  a  solution. 
They  indicate  that  each  customer  can  only  be  assigned  to  a  single  route  (2.3),  that 
all  vehicles  with  a  route  must  start  from  the  depot  (2.4),  that  all  edge  costs  are 
symmetrical  (2.5),  and  that  all  vehicles  must  return  to  the  depot  (2.6). 

Xijkif^ik  T  Si  T  tij  ijJ jk )  E  0  V k  G  AT,  (jl7  j)  G  A,  (2.7) 
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(2.8) 


Cj  ^  ^  %ijk  —  ^ ik  —  ^  ^  %ijk  ik  G  A ,  i'l  G  -k  - 

ie  A+(i)  jeA+p) 

ej  <  uj^  <  Vfc  G  A',  i  G  (0,  n  +  1),  (2.9) 

Equations  (2.7)-(2.9)  define  the  time  constraints  of  the  problem.  Equation  (2.7) 
indicates  that  the  arrival  time  at  location  i  plus  the  service  time  and  travel  time  to 
the  next  location  must  equal  the  arrival  time.  Equation  (2.8)  defines  the  need  for 
arrival  times  to  be  within  the  customers  time  window.  The  depot  also  has  a  time 
window  associated  with  it  (opening  and  closing  time)  which  all  vehicles  must  adhere 
to  (2.9). 


I\ 

T,d<  E  Xijk<^2ki  'ik  G  K,  (2.10) 

i£N  jeA+i  k= 0 

min  E  E  Cij^ijk  (2-11) 

kGK  (i,j)eA 

Equation  (2.10)  guarantees  that  the  vehicle  not  exceed  its  capacity  limit.  The 
objective  function  is  defined  by  Equation  (2.11)  which  illustrates  the  desire  to  min¬ 
imize  the  total  path  cost.  Note  that  the  objective  is  the  sum  of  the  path  costs  per 
vehicle  route.  The  path  cost  does  not  include  the  service  or  waiting  time.  Time  is 
only  a  constraint  that  causes  some  routes  to  be  infeasible. 

2. 3  Coalition  Forming 

Within  the  field  of  Artificial  Intelligence  and  Robotics  exists  a  field  of  study 
called  coalition  forming.  This  field  deals  with  the  problem  of  determining  the  optimal 
grouping  of  agents  to  different  targets  and  how  the  robots  should  get  there.  This  field 
of  study  is  a  precursor  to  the  development  of  the  Swarm  Routing  Problem.  In  Gerkey 
[20]  a  taxonomy  for  the  different  types  of  robot  tasking  problems  is  developed.  Gerkey 
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defines  the  tasking  problem  and  maps  it  to  a  combinatorics  set  covering  problem. 
By  doing  this  the  complexity  of  the  problem  is  defined  and  he  explores  to  define  a 
variety  of  solution  techniques.  The  problem  in  Gerkeys  research  termed  the  single- 
task  robots,  multi-robot  tasks,  time-extended  assignment  problem  is  the  same  premise 
that  is  proposed  in  this  research.  However  the  problem  is  not  expanded  by  mapping 
the  problem  model  from  the  VRPTW  as  is  done  here. 

2-4  Path  Planning 

The  process  of  path  planning  develops  a  sequence  of  steps  to  get  from  a  start¬ 
ing  point  to  an  end  point  based  on  the  internal  representation  of  the  terrain  to  be 
traversed.  The  path  is  then  optimized  for  a  minimal  cost.  Cost  to  get  to  a  point 
may  be  just  the  distance  traveled  but  can  include  fuel  used,  time,  or  hazards  within 
the  terrain.  A  path  planning  algorithm  is  a  method  used  to  calculate  a  path  plan 
given  knowledge  of  the  path  environment  and  a  set  of  conditions  or  constraints  that 
must  be  followed.  Once  generated,  the  plan  is  either  executed,  or  further  refined  to 
gain  improvement.  It  is  generally  assumed  that  once  the  plan  is  given,  the  machine 
becomes  autonomous  and  can  no  longer  interact  with  the  planner  [15]. 

Any  path-planning  algorithm  must  use  a  representation  of  the  terrain  in  order 
to  determine  what  the  cost  of  a  path  is.  Usually  this  is  a  grid  which  associates 
coordinates  with  locations  in  a  terrain  field.  The  points  of  interest  that  need  to 
be  visited  are  known  as  vertices  or  nodes.  Connectivity  between  these  vertices  is 
expressed  by  paths,  arcs  or  edges.  From  this  general  statement  of  the  purpose  of 
path  planning  the  concept  divides  into  two  realms  known  as  configuration  space  and 
trajectory  space. 

“Configuration  space”  is  a  problem  definition  where  all  possible  positions  of 
some  physical  system  are  defined.  The  solution  involves  determining  the  set  of  actions 
(torque,  translation,  etc.)  necessary  to  move  the  system  from  an  initial  state  to  a  goal 
state,  where  both  of  these  states  exist  within  the  defined  set  of  states  and  the  system 
never  goes  outside  the  set  of  possible  states.  Configuration  states  are  most  often  used 
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in  the  development  of  3D  robotic  motion  as  they  result  in  a  set  of  states  the  robot 
must  pass  through  to  reach  its  end  goal. 

“Trajectory  space”  also  deals  with  a  space  of  valid  and  invalid  positions  within 
a  space  of  possible  configurations.  The  result,  however  is  not  a  set  of  states  but  a  set 
of  trajectories  along  a  path  that  the  physical  system  passes  through  from  the  initial 
state  to  the  goal  state.  The  solution  is  determined  by  first  defining  a  path  between 
the  start  and  end  positions  that  passes  through  the  valid  and  invalid  positions.  The 
path  is  then  incrementally  adjusted  until  a  valid  or  optimum  path  is  formed.  These 
adjustments  are  made  by  determining  where  the  line  intersects  with  obstacles  and 
then  explores  the  subsection  of  the  path  to  determine  how  it  can  be  altered  to  avoid 
the  invalid  area. 

There  are  many  methods  for  dealing  with  both  of  these  problems  space  defini¬ 
tions.  For  small-scale  path  planning  problems  search  trees  are  often  used  to  determine 
the  optimum  path  the  physical  system  must  follow  [26] .  As  complexity  increases  into 
3-dinrensions  and  many  hazards,  evolutionary  methods  have  been  shown  to  be  effec¬ 
tive  for  determining  optimal  solutions  [25]. 

2. 4-0.1  Evolutionary  Solutions  in  Trajectory  Space.  One  of  the  first 
examples  of  evolutionary  computation  applied  to  trajectory  space  is  used  in  the  Evo¬ 
lutionary  Planner/Navigator  (EP/N)  created  by  Xiao  and  Michalewicz  [57].  The 
planner  treated  the  path  from  the  start  to  end  position  as  the  chromosome.  By  ma¬ 
nipulating  the  location  of  intermediate  points  through  the  genetic  operations  of  the 
Evolutionary  Algorithms  (EA)  the  path  was  evolved  into  an  optimal  path  consisting 
of  a  set  of  way-points.  This  work  was  later  refined  by  Sugihara  [46]  and  then  modified 
into  a  two  objective  problem  by  Castillo  [6].  These  previous  papers  defined  the  idea 
of  starting  with  a  straight  line  path  and  then  evolving  it  into  a  path  that  navigated 
around  obstacles  and  hazards  areas. 

In  order  to  deal  with  the  flight  dynamic  constraints  for  UAV  path  planning,  Slear 
used  a  concept  called  B-spline  curves  [44],  This  method  allowed  for  the  development 


18 


of  controllable  path  changes  while  the  path  is  being  evolved.  A  path  is  defined  by  sets 
of  tuples  (Xn,  Yn,  Zn)  where  each  set  consists  of  three  tuples.  By  modifying  which 
sets  are  placed  in  the  final  path  and  by  generating  random  tuple  sets  via  mutation 
optimal  paths  are  generated.  These  paths  can  be  formed  according  to  any  number  of 
objectives.  In  the  case  of  Slear  these  objectives  were  cost  and  risk  [44],  The  cost  of 
the  path  entails  the  fuel  required  to  climb  or  travel  the  distance  of  the  path.  The  risk 
is  the  proximity  to  enemy  locations  or  height  above  ground  (higher  altitudes  increase 
radar  visibility).  These  objectives  usually  do  not  complement  each  other,  as  a  low 
cost  path  generally  has  a  high  associated  risk  and  a  low-risk  path  generally  involves 
high  cost  (through  the  need  for  flying  low-over  uneven  terrain). 

The  domain  of  the  problem  consists  of  an  operational  n  x  m  grid  space,  a  terrain 
grid  G  (2.12)  and  a  location  set  L  (2.13).  The  location  set  is  the  set  of  targets  from 
the  VRPTW  described  in  Section  2.2. 


Ge(n-l)x(m-l)  (2.12) 

li  G  n  x  mV/  G  L  (2-13) 

The  problem  is  constrained  by  a  heading  change  requirement  defined  in  (2.14). 
This  constraint  indicates  that  changes  to  a  paths  heading  can  not  exceed  a  change 
greater  than  45°.  This  constraint  exists  to  ensure  that  solution  paths  do  not  require 
heading  changes  greater  than  a  UAV  is  capable  of  achieving.  This  constraint  is  of 
course  dependent  on  the  type  of  vehicle  being  modeled. 


Vpo--Rn  GP,  A6(pi,pi+i)  <  45° 

where  9  is  the  inbound  heading  at  pi 


(2.14) 
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The  objective  of  the  path  planning  problem  is  to  develop  a  minimum  cost  path 
P*  from  all  Z*  G  L  to  all  lj  G  L.  A  path  is  dehned  as  a  set  of  points  such  that 
pn  G  P.  The  determination  of  this  path  cost  is  found  in  five  parameters  divided  into 
two  separate  objectives  making  this  a  multi- objective  problem.  The  path  value  Qpath 
is  the  Euclidian  ground  distance  between  two  points.  The  climb  value  is  the 

change  in  altitude  that  occurs  over  some  path  length.  It  is  found  via  the  summation  of 
the  change  in  altitude  from  intermediate  points  in  the  path  being  examined.  Increases 
in  altitude  are  to  be  avoided  while  decreases  in  altitude  are  favorable  as  they  require 
less  power.  Detection  &  detect,  is  a  value  based  on  if  the  vehicle  is  in  a  known  detection 
area  and  for  how  long  which  can  also  be  thought  of  as  risk.  This  value  can  be 
calculated  as  either  a  binary  value,  which  would  indicate  that  detection  areas  are  to 
be  completely  avoided,  or  some  normalized  value  based  on  how  long  the  vehicle  is  in 
the  detection  area.  The  kill  cost  &kui  is  associated  with  being  in  range  of  the  enemies 
ability  to  fire.  This  is  similar  to  the  detection  area  but  is  generally  smaller  in  area. 
The  kill  cost  can  also  be  thought  of  as  a  binary  value  (i.e.  always  avoid)  or  some 
normalized  value  (i.e.  avoid  as  much  as  possible).  The  “terrain  objective”  terrain  is 
a  measure  of  how  many  points  the  UAV  can  be  seen  from.  Minimizing  this  objective 
ensures  the  UAV  can  not  be  seen  by  known  and  unknown  threats.  This  measure  is 
similar  to  the  altitude  cost  requirement,  however  by  inherently  desiring  low  altitude 
flight,  as  this  decreases  the  odds  of  being  spotted,  it  serves  to  effectively  contribute  to 
path  cost  optimization  even  though  it  is  calculated  differently.  These  five  objectives 
define  a  multi- objective  optimization  problem. 

2.5  Multi- objective  Evolutionary  Algorithms 

Evolutionary  and  stochastic  search  algorithms  have  been  applied  to  a  variety 
of  optimization  problems  with  great  success  [4,8,19,41],  A  review  of  the  concept  of 
evolutionary  algorithms  is  found  in  Appendix  A.  Evolutionary  algorithms  are  capable 
of  providing  polynomial  time  solutions  for  most  NP  and  NP-Complete  problems  [9] 
that  would  otherwise  require  an  exponential  or  intractable  solution  time. 
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NP  or  “Non-deterministic  Polynomial  time”  combinatorics  problems  [10]  are 
those  problems  where  finding  a  solution  requires  a  non-deterministic  Turing  machine 
and  is  intractable  otherwise.  The  solution,  once  found,  can  be  verified  in  polynomial 
time.  The  most  difficult  NP  problems  are  termed  NP-Hard,  which  refers  to  problems 
that  can  not  be  solved  through  exhaustive  search  in  either  tractable  or  intractable  time 
frames.  Some  NP-Hard  problems  can  be  transformed  into  NP-Complete  problems  if 
the  NP-Complete  solution  also  applies  to  the  NP-Hard  problem.  The  idea  is  that 
NP-Hard  problems  constitute  the  hardest  set  of  problems  to  solve  in  a  tractable  time 
due  to  an  unlimited  solution  space  while  NP-Complete  problems  comprise  a  subset  of 
those  problems  with  a  limited  solution  space  that  is  still  very  large.  Figure  2.2  shows 
the  overlapping  area  these  problem  types  exist  in.  A  classic  NP-Complete  problem 
is  the  TSP.  The  TSP  consists  of  a  set  of  n  points  on  a  grid  connected  by  weighted 
edges;  the  objective  is  to  determine  the  shortest  circuit  that  goes  through  all  n  points 
and  ends  where  it  began.  For  a  small  number  of  points  this  is  not  difficult,  however 
the  solution  complexity  time  of  the  problem  is  and  the  number  of  solutions  that 

must  be  analyzed  quickly  increases  to  an  intractable  number  [51]. 


NP-Hard 


NP-Hard 


P  =  NP 


P  =  NP 


Figure  2.2:  NP-Problem  subdivisions. 
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For  these  many  different  types  of  problems  a  variety  of  applicable  solution  tech¬ 
niques  exist.  Many  deterministic  solutions  exist  for  small  or  constrained  NP  problems 
which  rely  on  search  trees.  However  there  is  a  limit  to  how  effective  deterministic  (or 
even  heuristic  based)  search  methods  can  be,  especially  as  problem  size  and  complex¬ 
ity  grows. 

One  of  the  interesting  aspects  of  NP  problems  is  that  many  of  them  can  also 
be  modeled  as  multi-objective  problems.  There  are  two  conflicting  effects  of  a  multi¬ 
objective  problem.  The  first  is  that  the  problem  is  often  more  useful  because  it  more 
closely  approximates  reality,  but  this  comes  at  a  cost  of  the  second  effect:  increased 
complexity.  The  VRP  is  hardly  realistic,  however  the  VRPTW  described  earlier  in  this 
section  is  closer  to  reality,  and  if  more  constraints  were  applied,  such  as  heterogenous 
vehicles,  back-hauls  (pick  up  and  delivery),  or  multiple  depots,  the  problem  would 
become  even  more  realistic.  While  this  makes  the  solution  of  the  problem  more 
valuable  it  also  makes  an  optimal  solution  that  much  more  difficult  to  obtain. 

With  knowledge  of  the  effective  use  of  evolutionary  algorithms  and  the  need 
for  multi-objective  problem  solutions,  recent  work  has  focused  on  developing  multi¬ 
objective  evolutionary  algorithms  [3,33,62],  MOEAs  are  basically  the  same  as  a 
standard  GA  with  the  difference  of  how  solutions  are  evaluated  and  ranked.  In  a 
multi-objective  solution  the  concept  of  dominance  exists.  A  solution  is  said  to  dom¬ 
inate  if  there  is  no  other  solution  that  can  improve  one  of  the  objectives  without 
simultaneously  reducing  another  [9]. 

By  examining  the  dominance  of  different  solutions  and  ranking  them  accord¬ 
ingly,  a  more  accurate  selection  of  effective  solutions  can  be  made  for  future  gener¬ 
ations.  Also,  by  ranking  across  multiple  objectives  the  resulting  solution  achieves 
optimal  performance  across  all  objectives  without  being  biased  toward  any  one  ob¬ 
jective.  When  discussing  optimality  in  a  multidimensional  space  the  concept  of  the 
pareto  front  becomes  beneficial.  The  Pareto  front  is  the  set  of  non-dominated,  feasible 
solutions.  More  recent  MOEAs  [63]  [13]  use  this  understanding  of  the  Pareto  front  in 
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Figure  2.3:  Dominance  example  for  two  objectives  fl  and  f2  [56]. 

order  to  track  which  genotypes  are  developing  better  solutions.  This  concept  of  the 
Pareto  front  is  visualized  in  Figure  2.3.  Notice  that  the  resultant  set  of  solutions  in 
the  front  provides  solutions  with  different  tradeoff  values.  Which  solution  is  actually 
used  is  a  decision  made  by  a  user  or  by  some  pre-determined  rule. 

An  examination  of  the  literature  reveals  a  growing  appreciation  for  the  use 
of  MOEAs  in  complex  problems,  such  as  the  VRP  [48]  [29].  The  reason  for  this  is 
that  MOEAs  are  better  able  to  navigate  the  highly  irregular  solution  space  that  exists 
within  the  VRPTW.  What  has  also  been  shown  is  that  multi-objective  approaches  not 
only  develop  good  solutions  but  are  also  better  than  biased  single  objective  solutions 
for  the  optimization  of  any  of  a  problems  multiple  objectives  [33]  for  certain  problems. 

2.6  Self  Organized  UAV  Swarms 

The  idea  behind  self  organization  is  that  developing  the  control  structure  to 
make  a  group  of  UAVs  work  together  is  computationally  infeasible.  SO  theory  defines 
a  method  for  applying  simple  rules  to  individual  vehicles,  which  when  put  into  large 
groups  develop  the  emergent  behavior  relevant  to  the  mission  [7]. 

Within  UAV  simulation  a  control  methodology  must  be  established  that  defines 
how  each  UAV  flies  and  makes  decisions.  While  it  is  possible  to  manage  each  UAVs 
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Figure  2.4:  UAV  local  sensor  pattern  [37]. 

flight  with  a  deterministic  method  that  selects  optimal  flight  paths  this  is  an  extremely 
computationally  intensive  operation.  For  this  reason,  a  method  of  self-organized  con¬ 
trol  is  used.  Self-organization  [9]  refers  to  an  emergent  property  that  exists  within 
complex  systems  composed  of  very  simple  entities.  An  example  of  this  is  the  flocking 
of  birds.  Each  bird  has  no  global  knowledge  nor  understanding  of  direction,  they  only 
know  to  stay  within  a  certain  distance  of  its  neighbors  and  avoid  obstacles.  Figure  2.4 
illustrates  the  local  sensor  area  available  to  the  swarm  member.  The  net  effect  of  each 
bird  following  these  simple  rules  is  the  emergent  flocking  behavior  we  see  in  the  real 
world.  These  rules  were  defined  by  Reynolds  in  his  original  paper  on  the  subject  [38]. 

It  is  difficult  to  predict  what  low-level  rules  result  in  the  macroscopic  behavior 
desired  in  a  self-organized  entity.  For  this  reason  much  of  the  current  work  in  SO 
focuses  on  evolving  and  optimizing  existing  rule  sets  that  achieve  the  desired  behav¬ 
ior  [36] .  While  tuning  the  specific  parameters  is  difficult  and  time  consuming,  defining 
the  rules  is  not  as  complicated.  For  example,  in  a  swarm  of  UAVs,  the  application  of 
Reynold’s  three  original  flocking  rules  plus  a  rule  later  defined  in  [54]  for  flocking  can 
be  used.  These  rules  are  illustrated  in  Figures  2. 5-2. 8. 

The  first  two  items  can  actually  be  seen  as  a  single  rule  to  stay  within  a  certain 
range  of  nearby  flock  members.  Each  one  of  these  rules  represents  a  vector  of  influence 
for  the  swarm  member.  The  direction  the  swarm  member  flies  is  a  summation  of  these 
four  vectors.  The  result  of  these  four  rules  is  that  each  vehicle  avoid  crashing  into 
another,  all  members  fly  at  the  same  speed,  and  the  entire  flock  moves  toward  some 
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Figure  2.5:  Collision  avoidance:  avoid  hitting  nearby  swarm  members  [37]. 


Figure  2.6:  Swarm  centering:  stay  within  sight  of  nearby  swarm  members  [37]. 


Figure  2.7:  Velocity  matching:  match  velocity  of  nearby  swarm  members  [37]. 

objective  point.  The  basic  idea  of  self-  organization  is  not  that  complicated;  realistic 
simulations,  unfortunately,  are. 

There  are  quite  a  few  simulations  of  self  organized  behavior  available  online  [30] 
though  these  are  generally  limited  to  2D  simulations,  with  the  exception  of  the  Air 
Force  Institute  of  Technology  (AFIT)  UAV  simulator  and  a  few  others  [37].  The 
problem  with  conducting  simulations  is  the  inherent  complexity  of  managing  the 
global  information  of  the  swarm  and  its  environment  (0(n4)  for  the  4  basic  rules  per 
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Figure  2.8:  Migration:  maintain  velocity  to  a  set  way  point  [37]. 

UAV).  Recent  advances  in  computer  hardware,  data  structures,  and  parallel  software 
have  enabled  more  advanced  simulations  but  the  problem  complexity  is  still  very 
relevant,  especially  as  more  detailed  and  realistic  simulations  are  required. 

2. 7  Summary 

This  chapter  presents  a  short  history  of  previous  research  that  supports  this 
effort.  The  topics  of  the  VRPTW,  path  planning,  MOEAs,  and  SO  are  discussed  in 
their  relation  to  this  thesis  investigation.  Previous  research  is  also  explored  as  part 
of  topic  development.  The  problem  domain  consists  of  UAV  routing.  The  solution 
approaches  consist  of  MOEAs  (as  a  solution  to  the  routing  problem).  Next,  the  SRP 
is  defined  as  an  alternate  routing  model  and  along  with  the  VRPTW,  is  formulated 
as  a  multi-objective  optimization  problem. 
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III.  Problem  Definition 


This  chapter  formulates  the  SRP  in  a  manner  similar  to  that  of  the  VRPTW  in 
Chapter  II.  The  SRP  is  a  modification  of  the  VRPTW  that  allows  for  a  more 
realistic  routing  of  UAVs  to  targets.  The  chapter  concludes  with  a  multi-objective 
formulation  of  the  VRPTW  and  SRP. 

3.0.1  The  Swarm  Routing  Problem.  In  almost  all  variations  of  the  VRP,  or 
VRPTW,  it  is  assumed  that  all  vehicles  depart  from  the  depot  location  to  different 
targets  and  only  one  vehicle  visits  each  target  [51].  This  problem  model  is  appropriate 
because  each  vehicle  is  generally  assumed  to  be  ground  based.  However,  the  use  of 
a  UAV  swarm  introduces  an  interesting  aspect  which,  up  to  this  point,  has  not  been 
dealt  with.  When  dealing  with  a  swarm  of  UAVs  and  multiple  targets  to  visit,  it  is 
desirable  to  have  the  ability  to  route  the  swarm  between  targets  in  the  most  efficient 
manner  possible.  The  reason  for  this  is  that  often  many  targets  exist  on  a  battlefield 
which  need  to  be  visited  in  a  timely  manner,  while  also  utilizing  resources  in  the  most 
economical  fashion  possible.  It  would  be  much  more  efficient  to  be  able  to  consider 
UAVs  as  a  dynamic  group  rather  than  indivisible  units  that  can  only  visit  one  location 
at  a  time.  This  implies  an  imperative  to  take  advantage  of  the  divisibility  of  the  swarm 
and  route  subgroups  of  UAVs  to  different  targets,  as  it  is  deemed  efficient  to  do  so, 
and  then  regroup  at  other  targets. 

By  viewing  the  problem  in  this  manner  the  value  of  importance  changes  from 
edge  costs  between  targets  to  the  distance  traveled  by  each  UAV.  Within  the  VRPTW 
(and  VRPs  in  general)  each  vehicle  is  seen  to  have  some  capacity  associated  with  it 
that  is  used  to  satisfy  each  customer.  While  this  works  for  ground  based  delivery 
routing,  it  would  be  more  appropriate  to  view  target  satisfaction  as  the  number  of 
UAVs  on  target  in  some  time  window.  The  path  cost  associated  with  a  single  vehicle 
is  then  more  reflective  of  the  distance  that  needs  to  be  traveled,  and  the  use  of  many 
vehicles  capable  of  being  routed  through  multiple  targets  (causing  splits  and  joins  in 
the  group  in  the  process)  presents  a  more  realistic  and  useful  problem  model. 
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The  problem  construction  is  modeled  off  the  structure  of  the  VRPTW  found  in 
Toth  [51]  and  Tan  [47].  The  problem  domain  consists  of  a  network  G  =  (V,  A)  where 
V  =  {uoUi,  ■■■vn}  and  r>o  is  the  depot.  The  set  of  edges  is  defined  as,  A  =  {( Vi,Vj )  G 
V,i  7 -  j},  where  each  edge  has  associated  with  it  some  cost  c(vi,Vj).  The  cost  of  the 
edge  is  the  cost  of  travel  from  target  i  to  target  j.  For  now  we  assume  path  cost  as 
some  constant  travel  speed  for  each  UAV  making  the  travel  cost  simply  the  Euclidean 
distance  between  points. 

A  time  window  exists  for  all  customers  where,  E,  represents  the  earliest  start 
time  and,  L,  the  latest  arrival  time.  The  latest  arrival  time  is  the  point  at  which  the 
UAV  can  arrive  and  still  complete  the  service  time  defined  by  S.  If  the  vehicle  arrives 
earlier  than  E,  it  incurs  a  waiting  time,  W,  which  is  the  difference  between  the  arrival 
time  and  E.  The  total  time  a  vehicle  takes  to  complete  it’s  route  is  the  summation 
of  all  route  path  travel  costs,  waiting  times,  and  service  times  (Jj)  wi  +  s*)- 

The  total  path  time  must  not  exceed  the  latest  arrival  time  (i.e.  closing  time)  of  the 
depot. 

Each  customer  also  has  associated  with  it  some  demand,  D.  The  demand  is  an 
indication  of  the  number  of  UAVs  that  need  to  be  present  at  the  target  within  its  time 
window  for  the  required  service  time.  This  is  one  of  the  key  differences  between 
the  SRP  and  VRPTW.  Instead  of  demand  being  satisfied  by  a  UAVs  capacity, 
it  is  satisfied  by  the  number  of  UAVs  at  the  target.  The  service  time  indicates  the 
amount  of  time  the  UAVs  are  required  to  be  on  target. 

There  exists  K  homogenous  UAVs  each  of  which  has  some  travel  capacity,  F. 
and  unit  deliverable  capacity.  Groups  of  UAVs  are  classified  as  a  swarm.  A  swarm  can 
split  into  one  or  more  sub-swarms,  join  with  other  sub-swarms  into  a  larger  swarm, 
and  travel  along  path  edges  together  as  a  swarm.  It  is  assumed  that  join  and  split 
operations  only  occur  at  targets  in  order  to  simplify  problem  complexity.  The  travel 
capacity  is  not  a  deliverable  value  as  it  has  been  in  previous  versions  of  this  problem, 
it  is  only  an  indication  of  how  far  the  individual  UAV  can  fly.  This  constraint  can 


be  viewed  as  the  UAVs  power  supply  limitation.  The  deliverable  value  that  satisfies 
the  target  is  equal  to  the  total  number  of  UAVs  present  in  a  given  location  at  a  given 
time.  This  value  fulfills  the  demand  requirement  of  the  target  during  its  service  time. 

The  solution  to  the  problem  is  the  same  as  the  VRPTW,  a  list  of  ordered  targets 
for  each  vehicle,  such  that  the  visitation  to  each  target  fulfills  all  target  needs  without 
violating  any  time  or  demand  constraints.  Note,  that  the  cost  of  a  route  is  not  the 
total  time  the  route  takes  to  complete,  it  is  only  the  sum  cost  of  the  edges  the  vehicle 
traverses.  The  objective  remains  the  same:  determine  the  set  of  paths  for  the  UAVs 
such  that  the  total  distance  is  minimized. 

The  following  is  a  mathematical  programming  formulation  of  the  SRP,  based 
on  the  definition  found  in  Toth  [51].  Vehicles  are  defined  within  the  problem  by  their 
inclusion  in  a  flow  variable,  which  is  a  binary  value  indicating  if  UAV,  k,  exists 
on  the  path  that  connects  (i,j)  G  V  at  any  point  in  the  solution.  A  time  variable,  a 
indicates  the  start  time  of  UAV  k  at  location  i.  The  subscript  j  G  A^z)  indicates  the 
set  of  edges  from  i  to  j  where  j  is  not  equal  to  i,  the  plus  or  minus  indicates  either  a 
forward  or  backward  move  along  the  path. 

Aij  -  Edge  cost  between  i  and  j 

Vn  -  Network  vertices  for  n  target  (no  is  the  depot) 

En  -  Earliest  arrival  time  of  target  n 

Ln  -  Latest  arrival  time  of  target  n 

Sn  -  Service  time  of  target  n 

Dn  -  Demand  of  target  n 

K  -  Set  of  UAVs 

Efc  -  Travel  capacity  of  UAV  k 


%ijk  G  {0,  1} 


\/k  G  K,  (i,j)  G  A,  (3.1) 
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^ ik  ^  0 


i  G  N,k  G  K, 


(3.2) 


Equations  (2.1)  and  (2.2)  define  the  flow  and  time  variables  used.  The  flow 
variable  is  a  binary  value  that  indicates  vehicle,  k,  travels  from  location,  i  to  j, 
if  equal  to  one,  and  zero  otherwise.  The  time  variable  specifies  the  start  time  at 
location,  i,  by  vehicle  k. 


y  Xijk  =  1  Vi  £  iV,V/c  £  K,  (3.3) 

jeA+p) 

y  X0jk  =  1  Wk  £  K,  (3.4) 

j'eA+(o) 

y  xi:n+  i)fe  =  1  V/c  £  K,  (3.5) 

ieA~(n+l) 

y  Xijk-  y  Xijk  =  0  \/k  e  K,\/i  e  N,  (3.6) 

ieA+(j)  ieA-(j) 


Equations  (3.3)-(3.6)  define  the  edge  constraints  of  the  graph  in  a  solution. 
They  indicate  that  the  each  vehicle  visits  a  customer  only  once  (3.3),  that  all  vehicles 
must  start  from  the  depot  (3.4),  that  all  edge  costs  are  symmetrical  (3.5),  and  that 
all  vehicles  must  return  to  the  depot  (3.6).  Note,  in  Equation  (3.3)  the  lack  of  the 
summation  over  K  which  removes  the  constraint  that  each  customer  be  visited  by 
only  a  single  vehicle. 


Xijk{i^ik  ~\~  Si  T  tjj  uj jk^)  0  V/c  £  K)  (z,  j)  £  zl,  (3.7) 


Si  ^  (  Xijk  ^  uj ik  ^  li  ^  (  Xijk  V/c  £  ./V,  Vz  £  TV,  (3.8) 

jeA+p)  jeA+p) 
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C-i  —  ^ ik  —  l i 


Wk  G  K ,  %  G  (0,  n  +  1), 


(3.9) 


Equations  (3.7)-(3.9)  define  the  time  constraints  of  the  problem.  Equation  (3.7) 
indicates  that  the  arrival  time  at  location  i  plus  the  service  time  and  travel  time  to  the 
next  location  must  equal  the  arrival  time  at  the  next  customer.  Equation  (3.8)  defines 
the  need  for  arrival  times  to  be  within  the  customers  time  window.  The  depot  also 
has  a  time  window  associated  with  it  (opening  and  closing  time)  which  all  vehicles 
must  adhere  to  (3.9). 


K 

Xjjk  >  T!  h  Vfc  G  K,  (3.10) 

i£N  jeA+i  k= 0 

E  CijXijk  <  Fk  Wk  e  K,  (3.11) 

jeA+i 

Up  to  this  point  the  formulation  is  basically  been  the  same  as  the  VRPTW  with 
the  exception  of  Equation  (3.3).  Equations  (3.10)-(3.11)  are  what  separate  the  SRP 
from  the  VRPTW.  Equation  (3.10)  indicates  that  the  demand  of  each  customer  is 
satisfied  by  the  number  of  vehicles  on  location  and  that  that  number  must  be  either 
equal  to  or  greater  than  the  demand  of  the  target.  Equation  (3.11)  indicates  that  the 
total  cost  of  the  path  for  a  vehicle  not  exceed  the  vehicles  flight  cost  limit. 


min  E  E  CijXijk  (3.12) 

k£K  ( i,j)£A 

The  single  objective  function  is  defined  by  Equation  (3.12)  which  illustrates  the 
desire  to  minimize  the  total  path  cost  for  all  vehicles.  The  path  cost  does  not  include 
the  service  or  waiting  times.  Time  is  only  a  constraint  that  causes  some  routes  to  be 
infeasible,  the  cost  is  the  total  distance  traveled. 
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This  formulation  has  introduced  the  SRP  as  a  modification  of  the  VRPTW. 
By  changing  the  constraints  of  customer  visitation  and  how  a  customers  demand  is 
satisfied  the  problem  becomes  a  more  realistic  model  for  routing  UAVs  to  multiple 
targets  within  a  time  window.  Up  to  this  point  problem  and  the  VRPTW  defined 
in  Chapter  II  have  only  been  single  objective  formulations.  In  the  next  section  the 
problems  are  reviewed  in  terms  of  multiple  objectives. 

3.1  Multi- objective  Formulation  for  VRPTW  and  SRP 

In  Section  2.2  the  VRPTW  is  defined  and  in  the  previous  section  a  variant,  the 
SRP,  is  defined.  The  objective  functions  for  these  two  problems  indicates  that  only 
path  length  is  of  critical  interest.  Even  though  path  length  is  a  primary  objective 
it  is  not  the  only  objective  that  can  be  optimized  for  in  the  solution.  Consider  the 
following  situations  that  may  occur  within  a  problem: 

•  Vehicle  exceeds  its  capacity  to  serve  a  route  -  when  this  happens  the 
route  can  be  split  into  2  or  more  routes.  The  split  is  made  when  the  customer 
demand  causes  the  capacity  of  the  vehicle  to  be  exceeded.  If  the  original  route 
(1,2, 3, 4, 5, 6)  causes  the  vehicle  capacity  to  be  exceeded  at  customer  4,  the  route 
is  split  into  two  routes  (1,2,3)  and  (4,5,6)  and  a  new  vehicle  is  added  to  the 
solution.  This  would  be  repeated  in  all  routes  until  the  capacity  is  met. 

•  Vehicle  arrives  early  to  a  customer  -  when  this  happens  the  service  time 
for  the  customer  is  increased  by  the  time  spent  waiting 

•  Vehicle  violates  a  time  window  by  arriving  late  -  a  new  vehicle  and 
route  are  added,  splitting  the  route  as  described  when  capacity  constraints  are 
encountered. 

From  these  situations  we  see  that  the  solution  to  alleviating  capacity  and  time 
violations  is  to  increase  the  number  of  vehicles  (and  routes).  Increasing  the  number 
of  vehicles  is,  of  course,  regressive  to  the  development  of  optimal  paths  lengths.  This 
is  due  to  the  introduction  of  depot  travel  times  for  each  new  route.  Every  time  a  new 
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route  is  added  the  vehicle  must  first  travel  from  the  depot  to  a  location,  making  the 
addition  of  new  vehicle  routes  generally  cause  an  increase  in  total  path  length  (though 
not  always).  It  is  therefore,  advantageous  to  define  3  objectives  to  minimize  for:  path 
length,  vehicle  count  and  total  wait  time.  By  optimizing  for  these  three  objectives 
we  seek  solutions  with  complementary  aspects,  and  better  solutions  overall.  The 
objective  functions  now  consist  of  the  following  equations: 

min  E  E  CijXijk  (3.13) 

k£K  (■ i,j)eA 

I< 

min  k  (3-14) 

k= 0 

min  E  E  tik^ijk  (3.15) 

keK  ( i,j)eA 

Equation  (3.13)  defines  the  minimization  of  the  path  length.  Equation  (3.14) 
defines  the  minimization  of  the  number  of  vehicles  used.  Equation  (3.15)  defines  the 
minimization  of  the  wait  time  where  the  variable  t  is  defined  in  equation  (3.16)  for 
the  VRPTW  and  (3.17)  for  the  SRP.  The  wait  time  for  the  SRP  is  defined  as  the 
difference  between  a  vehicles  arrival  time  and  the  arrival  time  of  the  latest  vehicle,  if 
the  latest  arriving  vehicles  time  is  past  the  earliest  arrival  time  of  the  customer.  The 
latest  arriving  time  is  used  because  operation  on  a  customer  can  not  begin  until  all 
vehicles  are  present.  These  objective  functions  apply  to  both  the  VRPTW  and  SRP. 

{e-i  -  wik,  if  E  >  wik 
0,  otherwise 
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tik  < 


et  -  wik,  if  E  >  wik  n  E  >  wiu 
wiu  -  wik,  if  wiu  >  E 
0,  otherwise 

where  u  is  the  latest  arriving  vehicle 


(3.17) 


3.2  Purpose  in  multi- objective  formulation 

Optimizing  across  multiples  objectives  is  done  with  two  purposes  in  mind.  First, 
to  support  the  idea  that  a  multi-objective  formulation  is  capable  of  navigating  the 
solution  space  more  effectively  than  optimizing  for  only  a  single  objective.  Since  the 
objectives  complement  each  other  it  would  seem  logical  that  optimizing  over  all  of 
them  would  achieve  better  results.  What  is  also  noteworthy  is  that  by  optimizing  for 
these  different  objectives,  solutions  with  decreased  path  lengths  should  be  found  as 
opposed  to  optimizing  solely  for  path  length.  This  idea  was  first  proposed  and  tested 
in  Ombuki  [33]  with  beneficial  results  (i.e.  benchmark  values  were  not  made  worse 
from  the  multi-objective  approach  compare  to  the  single  objective  approach). 

The  reason  multi-objective  formulation  is  more  effective  is  because  the  problem 
under  consideration  has  such  an  irregular  solution  space.  Time  constraints  introduce 
irregularities  to  the  Pareto  front  such  that  non-dominated  solutions  become  more 
isolated.  This  problem  is  exacerbated  as  more  constraints  are  applied  to  the  prob¬ 
lem  such  as  heterogenous  vehicle  fleets  or  pick-up  and  delivery  problems.  Only  by 
optimizing  across  multiple  objectives  can  the  solution  space  be  traversed  accurately 
enough  to  allow  the  determination  of  non-dominated  solutions.  We  now  proceed  with 
the  high  level  solution  design  to  this  fully  formed  multi-objective  problem. 


3. 3  Summary 

This  chapter  introduced  the  SRP  routing  problem  variation  as  a  model  for 
routing  swarms  of  UAVs.  By  eliminating  the  restriction  of  single  city  visitation  for 
each  vehicle  and  changing  target  satisfaction  to  be  based  on  vehicle  count  the  SRP 


34 


serves  as  a  better  problem  model  for  UAV  routing.  A  multi-objective  formulation 
for  the  both  the  VRPTW  and  SRP  shows  the  objectives  can  be  optimized  for:  path 
length,  vehicle  count,  and  wait  time  per  vehicle.  The  next  chapter  introduces  the 
high  level  solution  design  for  these  two  problems. 
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IV.  High  Level  Concept  Design 


This  chapter  presents  the  high  level  design  that  forms  a  solution  to  the  multi¬ 
objective  VRPTW  and  SRP  defined  in  Chapter  III.  An  analysis  of  problem 
complexity  precedes  a  development  of  the  algorithmic  solution.  This  development 
follows  a  process  of  determining  the  type  of  solution,  the  structure  of  the  algorithm 
itself,  and  the  individual  components  that  make  up  its  design. 

4-1  Design  Objectives 

The  design  objectives  for  this  research  are  to  develop  a  solution  procedure  for 
the  multi-objective  routing  problems  described  in  the  previous  chapter.  This  solution 
must  return  information  that  can  be  integrated  with  previous  work  on  UAV  path 
planning  [44]  and  simulation  [11],  The  end  result  of  this  design  is  a  fully  developed 
form  of  an  algorithm  to  be  applied  to  the  VRPTW  and  SRP  whose  output  is  a  set  of 
feasible  vehicle  routes. 

A  discussion  of  the  problem  complexity  is  performed  in  order  to  illustrate  design 
requirements  of  the  algorithm.  An  introduction  to  the  form  of  the  required  solution 
is  then  shown;  this  form  requires  two  aspects:  a  selection  method  and  set  of  genetic 
operators.  A  description  of  two  multi-objective  EAs  is  given  in  support  of  the  first 
requirement.  This  is  followed  by  a  set  of  genetic  operators  that  compose  the  GA 
solution  for  the  VRPTW  and  SRP,  each  contained  in  its  own  respective  section  1 . 

4-2  Problem  Complexity 

The  VRPTW  is  NP-Hard  [10,  51]  meaning  that  no  polynomial  time  solution 
exists  unless  V  =  AfV.  This  complexity  valuation  is  determined  by  analysis  of  the 
VRP  as  NP-Hard  [45]  and  through  restriction,  determining  that  the  VRPTW  is  also 

lrThe  term  evolutionary  algorithm  and  genetic  algorithm  are  almost  interchangeable  and  often 
used  as  such.  Evolutionary  algorithm  is  usually  considered  a  broader  term  encompassing  many 
different  type  of  evolutionary  computation.  Genetic  algorithms  refer  exclusively  to  optimization 
and  search  algorithms  which  employ  chromosomal  solution  manipulation,  where  the  chromosome  is 
directly  related  to  the  solution  (this  excludes  genetic  programming) 
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NP-Hard.  For  a  fixed  number  of  vehicles  the  problem  becomes  NP-Complete  [51]  for 
either  problem  variation.  The  SRP,  as  an  extension  of  the  VRPTW,  can  likewise  be 
classified  as  NP-Hard  in  the  general  sense  and  NP-Complete  for  fixed  fleet  sizes. 

To  determine  the  solution  space  complexity  we  start  with  the  complexity  of  a 
restrained  problem,  called  the  TSP.  In  TSP  the  number  of  vehicles  is  constrained 
to  one,  making  the  size  of  the  problem  1/2 (n  —  1)!  for  n  >  2,  where  n  is  the  num¬ 
ber  of  customers/targets  in  the  problem.  The  solution  space  size  of  the  VRP,  S,  is 
approximated  by  Equation  (4.1). 


S 


8n\/3 


1)! 


(4.1) 


This  equation  is  found  by  combining  an  integer  partition  distribution  generating 
function  and  the  the  complexity  of  a  single  n  customer  routing  problem  [55].  The 
result  of  this  approximation  yields  a  solution  space  complexity  of  (P(expn  n\).  For 
the  VRPTW  the  solution  space  is  the  same  as  it  contains  the  same  number  of  total 
solutions  as  the  VRP  minus  some  constant  number  of  solutions  made  invalid  by  time 
constraints.  This  same  idea  would  apply  to  any  VRP  variation  (i.e.  the  SRP  solution 
space  complexity  is  that  of  the  VRPTW).  This  rapid  increase  in  the  solution  space  is 
the  source  for  the  VRP  being  such  a  difficult  problem  to  solve  in  the  general  sense. 

The  complexity  of  the  VRPTW  leads  to  the  conclusion  that  a  polynomial  time 
solution  most  likely  does  not  exist.  The  best  course  of  action  is  to  then  pursue  a 
heuristic  based  solution  most  likely  structured  within  a  stochastic  algorithm  [9].  This 
assumption  is  backed  by  published  research  on  the  success  of  stochastic  solutions  to 
the  VRP  and  VRPTW  [21]  [34]  [47]  [3], 


4-3  Multi-objective  Algorithm  Design  Choices 

Many  applicable  algorithm  choices  exist  for  solving  the  VRPTW  and  most  are 
capable  of  returning  optimal  solutions  in  a  reasonable  time  as  seen  by  the  choices 


37 


examined  in  this  section.  The  choices  have  been  constrained  to  stochastic  based 
algorithms  as  our  previous  complexity  analysis  has  shown  that  a  deterministic  search 
process  would  be  intractable. 

4-3.1  Ant  Colony  Optimization.  Ant  colony  optimization  is  defined  by  a  set 
of  artificial  ants,  or  agents,  navigating  through  a  defined  search  space  attempting  to 
optimize  some  problem  solution.  The  route  the  ants  construct  within  the  search  space 
forms  the  solution.  Each  ant  exists  at  some  point  in  the  solution  space  and  makes  a 
local  decision  of  where  to  go  next.  This  decision  is  inherently  stochastic  but  is  based 
on  the  idea  of  the  pheromone  matrix.  In  nature,  ants  lay  pheromone  as  they  wander 
through  a  space.  These  pheromones  attract  other  ants  to  follow  the  same  path.  As  the 
ants  encounter  desirable  factors  more  ants  start  to  follow  the  same  path  reenforcing 
the  path.  This  idea  is  presented  visually  in  Figure  4.1.  While  this  does  not  lead  to 
an  optimum  solution  (true  best  solution)  it  does  often  lead  to  optimal  solutions  (high 
value  that  may  or  may  not  be  the  true  best  value),  especially  in  Euclidian  spaces 
(specific  to  routing  problems). 


Figure  4.1:  ACO  Paths  Being  Created. 

The  decision  of  which  direction  an  ant  travels,  within  the  ACO  model,  is  defined 
by  the  probability  definition  in  Equation  (4.2).  The  equation  is  the  pheromone  inten- 
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sity,  r,  times  the  distance  between  two  points,  divided  by  the  summation  of  all  other 
options  calculated  the  same  way.  77  is  a  visibility  term  that  denotes  the  inverse  of 
the  distance  between  two  points.  This  causes  the  ant  to  have  a  preference  for  shorter 
routes.  This  probability  is  a  function  of  the  pheromone  matrix  which  is  updated  as 
solutions  are  found.  Altering  the  update  method,  seen  in  Equation  (4.4),  to  update 
under  certain  conditions  or  only  for  the  highest  performing  ant  also  impacts  the  per¬ 
formance  of  the  algorithm.  The  pheromone  matrix  is  also  altered  by  Equation  (4.3) 
where  p  is  the  evaporation  rate  [0  —  1] .  The  net  result  of  these  methods  is  to  enforce 
the  use  of  paths  that  have  lead  to  good  solutions.  Algorithm  1  incorporates  these 
ideas  into  a  pseudo  code  representation. 


Algorithm  1  ACO  Algorithm 
1:  procedure  ACO([6esf]  =  max Jt, N, To) 

2:  initialize  rtj  //  usually  initialized  with  the  same  To 

3:  initialize  best 

4:  place  each  ant  fcona  randomly  selected  edge 

5:  t  <-  1 

6:  while  t  <  max  At  do 

7:  for  i  =  1  to  N  do 

8:  Build  a  solution  by  applying  a  probabilistic 

9:  transition  rule  (e  —  1)  times. 

10:  The  rule  is  a  function  of  r  and  r) 

11:  e  is  the  number  of  edges  on  the  graph  G 

12:  end  for 

13:  Evaluate  each  solution 

14:  if  an  improved  solution  is  found  then 

15:  update  the  best  solution  found 

16:  end  if 

17:  Update  pheromone  trails 

18:  t  < —  t  T  1 

19:  end  while 

20:  end  procedure 
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(4.2) 


V  (hj)  =  Vd{i,j),P  =  2 


Tij  <-  (1  -  p)Tij  +  A Tij 


(4.3) 


Supd  arg  max  F(s) 


(4.4) 


This  is  an  example  of  a  popular  elitist  update  strategy  where  for  every  iteration 
in  the  number  of  generated  solutions  the  one  with  the  highest  fitness  value  is  used  to 
update  the  matrix. 


4- 3. 1.1  Applying  ACO  to  the  MOP  VRPTW.  ACO  can  be  applied 


to  multi-objective  problems  in  one  of  two  ways.  A  set  of  different  colonies  can  be 
run  simultaneously  attempting  to  optimize  each  objective  before  forming  a  final  so¬ 
lution  [19],  or  a  single  colony  can  be  run  where  the  probability  equation  described 
previously  (Equation  4.2)  is  based  not  only  on  the  pheromone  matrix  but  also  on  ob¬ 
jective  weights  that  become  associated  with  the  pheromone  matrix  [3].  Each  option 
has  advantages  and  disadvantages.  Attempting  to  solve  for  all  objectives  simultane¬ 
ously  has  the  risk  of  creating  bias  toward  some  objectives,  and  using  different  colonies 
has  a  higher  computational  complexity. 

Optimizing  for  all  objectives  simultaneously  requires  each  of  the  objectives  de¬ 
scribed  in  Section  3.1  to  be  incorporated  into  the  fitness  evaluation  used  to  calculate 
the  quality  of  a  solution.  Baran  suggests  a  method  where  each  ant  initially  creates 
a  feasible  solution  by  adding  cities  to  routes  as  they  are  deemed  feasible  (within  the 
time  window  and  vehicle  capacity)  [3].  From  this  point,  a  variation  of  the  pheromone 
matrix  is  used  to  incorporate  the  value  of  pheromone  paths  to  the  different  objectives 
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being  optimized  (4.5).  This  new  equation  introduces  the  r\j  value  which  indicates 
time  between  cities.  In  this  way,  a  single  colony  can  be  used  to  develop  a  solution 
that  is  optimal  for  all  objectives. 


Pr  (i,j) 


T.allawed^ 

1  /d(i,j),rjj(i,j)  =  1  /t(i,j),@  =  2 


(4.5) 


4-3.2  Particle  Swarm  Optimization.  Particle  Swarms  are  multi  agent  sys¬ 
tems  inspired  by  biological  systems.  They  are  a  progression  in  the  area  of  natural 
computing  borrowing  ideas  from  genetic  algorithms,  ant  colonies,  and  self  organiza¬ 
tion.  The  general  algorithm  consists  of  a  population  of  agents  in  D-dimensional,  real 
valued  space  (9T°).  These  agents  move  through  the  space  attempting  to  locate  opti¬ 
mal  solutions  (as  seen  in  Figure  4.2).  Where  a  given  agents  movement  is  based  on  its 
own  experience,  its  current  location,  and  the  experience  of  the  agents  around  it. 

The  agents  in  the  swarm  are  defined  as  Xt  =  (. xu ,  X&,  Xu, ...,  xto),  where  X  is  a 
position  in  the  space.  Each  agent  also  has  an  associated  velocity  Vi  =  (va,  Vn, ...,  vw). 
The  velocity  a  particle  travels  is  defined  by  Equation  (4.6).  The  variables  Cj  represent 
positive  constants  while  ip  represents  a  uniform  random  number  generation  function. 


Figure  4.2:  Particles  on  a  Solution  Space. 
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The  variable  w  represents  the  inertia  of  the  problem  which  corresponds  to  the  expe¬ 
rience  the  agent  has  had  in  the  search  space.  The  p  term  is  the  best  value  found  so 
far  and  the  best  found  by  a  neighbor  (pgd)-  These  terms  work  to  move  the  particles 
velocity  toward  better  solutions.  These  best  values  can  be  set  to  either  a  local  best 
or  a  global  best.  The  impact  of  either  depends  on  the  problem  being  solved,  though 
global  knowledge  is  usually  more  beneficial  [58]. 


Vid(t  +  1)  =w  x  vid(t)+  (4.6) 

Cl  X  pi  X  \pid(t)  -xid^t)(t)}  + 
c2  x  p2  x  [pgdit)  -  Xid(t)(t)] 

Each  particle  also  has  a  max  and  min  velocity  that  it  can  achieve.  This  value  is 
set  by  the  user  and  keeps  the  particle  from  either  shutting  down  or  accelerating  out 
of  the  solution  space.  These  ideas  are  incorporated  into  pseudo  code  in  Algorithm  2. 

4-3.2. 1  Applying  PSO  to  the  MOP  VRPTW.  In  Zhu  [61]  and  Yan- 
wei  [58],  a  Particle  Swarm  Optimization  (PSO)  algorithm  is  applied  to  a  VRPTW.  A 
solution  is  structured  asa2xn  +  fc  -  1  particle,  where  n  is  the  number  of  customers 
visited  by  k  vehicles.  Thus  a  solution  to  a  six  dimension  problem  (with  zero  as  the 
depot)  would  be: 


Node  :00123456 
Index  :  26874315 

The  objectives  to  be  minimized  are  the  total  path  length,  the  number  of  vehicles 
used,  and  the  time  lost  waiting  for  customers.  In  Zhu  [61],  infeasible  solutions  (those 
that  violate  capacity  or  time  constraints)  are  penalized  instead  of  being  repaired  or 
thrown  out.  Thus,  each  objective  constraint  has  some  weight  associated  with  it. 


42 


Algorithm  2  PSO  Algorithm  from  [7] 


1 

2 

3 

4 

5 

6 
7 


procedure  [A"]  =  PSO(maxJt,ACi,AC'2,  vmax,vmin) 
initialize  A"  //  usually  a:,,  Vi,  is  initialized  at  random 
initialize  A  X;  //  at  random,  A  X;  G  [t’mm,umax] 
t  <—  i 

while  t  <  max  jit  do 
for  i  =  1  to  A  do 

if  0(x  i)  <  g(pi)  then 


8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
23 


Pi  =  Xi 

end  if 

g  =  i  / /  arbitrary 

//  for  all  neighbors 

for  j  =  indexes  of  neighbors  do 

if  9(Pj )  >  g(pg) then 

g  =  j  //  index  of  best  neighbor 

end  if 
end  for 

A  Xi  ^A  Xi  +  p  ®  (pi  -  Xi)  +  ip  <S>  ( pg  -  Xi) 
A  Xi  G  [Umin,  ®moi] 

Xi  <—  Xj+  A  Xi 

end  for 

t<-t  +  l 

end  while 
end  procedure 


These  values  are  incorporated  into  the  fitness  function  that  determines  how  good  a 
solution  is  and  ultimately  where  a  particle  goes. 

Methods  for  incorporating  multi-objective  capabilities  include  running  multiple 
populations  for  individual  objectives  and  then  combining  the  best  results  for  a  final 
search  [3]  and  using  a  memory  system  to  store  non-dominated  solutions  [16].  In 
Hu  [22],  a  method  of  randomly  selecting  non-dominated  points  to  be  incorporated 
into  the  population  is  used  to  ensure  exploration  of  the  search  space  as  well  as  pursuit 
of  all  objectives.  The  implementation  of  such  an  algorithm  would  closely  resemble 
the  general  implementation  in  Section  4.3.2.  Since  PSO  algorithms  do  not  rely  on 
structure  like  GAs  the  only  aspect  that  is  problem  specific  is  the  solution  structure 
and  the  fitness  function.  The  search  of  the  solution  space  is  carried  out  by  a  vector 
definition  in  Equation  (4.6),  defined  by  fitness  values,  not  the  solution  structure. 
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4.3.3  Evolutionary  Computation.  For  a  review  of  basic  evolutionary  compu¬ 
tation  concepts  please  refer  to  Appendix  A.  Evolutionary  computation  works  on  the 
idea  that  existing  problem  solutions  can  be  modified  in  an  intelligent  manner  in  order 
to  allow  for  the  creation  of  new  and  better  solutions.  Evolutionary  Computation  (EC) 
methods  are  very  applicable  to  routing  problems  because  routing  solutions  can  be  cod¬ 
ified  in  very  simplistic  chromosome  structures.  This  simplification  allows  for  easier 
construction  of  the  algorithm  and  a  faster  run  time.  The  next  section  presents  an 
expanded  description  of  the  multi-objective  GA. 

4-4  Multi- objective  Genetic  Algorithm  Design 

Multi-objective  GAs  differ  from  single  objective  GAs  in  how  solutions  are  strat¬ 
ified.  In  a  single  objective  algorithm,  determining  the  quality  of  a  solution  is  a  simple 
matter.  If  the  objective  is  minimization  it  is  a  simple  compare  operation  to  the  low¬ 
est  value.  For  multi-objective  optimization  this  process  is  not  as  simple.  While  it 
is  possible  to  weight  objectives  in  order  to  obtain  a  single  value  associated  with  the 
solution  this  is  not  an  advisable  pursuit.  Weighted  objectives  introduce  bias  because 
no  weighting  procedure  can  accurately  treat  the  different  objectives  in  a  manner  such 
that  all  objectives  are  optimized  effectively  [9].  To  accurately  classify  solutions  over 
a  multi-objective  domain  ,a  ranking  procedure  must  be  used  that  takes  into  account 
not  only  the  value  of  the  solution  across  the  different  objectives  but  also  its  proximity 
to  other  solutions,  which  is  an  indication  of  the  value  of  the  information  the  solution 
contains.  Chapter  II  expands  on  this  idea  of  pareto  front  searching. 

There  are  three  major  components  to  the  EA  design:  the  replacement  method, 
the  chromosome  structure,  and  the  genetic  operators.  The  replacement  method  de¬ 
termines  which  solutions  are  kept  after  a  new  generation  is  created.  Within  this  step 
the  solutions  are  ranked  and  discarded.  The  selection  methods  chosen  are  shown  in 
the  context  of  the  complete  EA.  The  sections  following  the  replacement  method  ex¬ 
pand  on  the  chromosome  structure,  population  initiation,  and  genetic  operators.  The 
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chromosome  structure  is  a  critical  step  which  drives  the  effectiveness  of  the  entire 
algorithm  and  how  the  different  genetic  operations  function. 

4-5  Replacement  Method 

In  this  section  two  replacement  methods  are  shown  in  the  context  of  the  GA 
they  form.  Both  of  these  methods  rely  on  non-dominated  sorting  and  objective  space 
distance  to  rank  and  select  solutions  for  the  next  generation.  How  Pareto  ranking  and 
dominance  is  utilized  multi-objective  search  is  discussed  in  Section  2.5.  As  arguments 
can  be  made  for  any  given  selection  method  for  any  given  problem  two  algorithms 
were  selected  so  that  their  results  could  be  statistically  compared.  While  many  dif¬ 
ferent  MOEA  selection  methods  exist,  Non-dominating  Sorting  Genetic  Algorithm 
2  (NSGA2)  and  Strength  Pareto  Evolutionary  Algorithm  2  (SPEA2)  were  chosen  for 
their  general  acceptance  within  the  research  community  and  previous  experiments 
showing  their  effectiveness  [35]  [62], 

4-5.1  Non- dominating  Sorting  Genetic  Algorithm  II.  NSGA2  uses  an  elitist 
sorting  mechanism  of  the  non-dominated  points  to  first  organize  the  solution  set  [14]. 
The  result  of  this  sorting  is  a  set  of  solution  ranks.  The  first  rank  is  the  hard  non- 
dominated  set.  Hard  non-dominated  refers  to  a  point  that  is  not  dominated  by  any 
other  solution,  as  apposed  to  soft  non-dominated  solutions  which  are  only  dominated 
by  those  points  in  the  first  rank.  Each  decreasing  rank  is  dominated  by  more  points. 
These  points  are  then  compared  to  each  other  in  order  to  determine  the  distribution  of 
points  in  the  current  Pareto  front  and  which  points  contribute  best  to  an  exploration 
of  the  solution  space.  This  process  is  called  crowding-distance-assignment  and  is 
shown  in  Algorithm  3. 

The  notation  Y[i\:m  refers  to  the  m-th  objective  value  for  the  the  f-th  solution. 
By  using  the  crowding  distance  and  ranking  procedure  the  solutions  are  ordered  by 
how  ’’good”  they  are.  A  visualization  of  this  idea  is  shown  in  Figure  4.3.  The  next 
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Algorithm  3  Crowding  Distance  Assignment 
l:  procedure  CrowdDis Assign  (T) 

2:  1=  |r| 

3:  for  each  i,  T[i\distance  =  0  do 

4:  for  each  objective  m  do 

5:  r  =  Sort(r,m) 

6:  T\l]  distance  =  T\l]  distance  =  oo 

7:  for  i  -  2 toil  -  1)  do 

8:  T[i\distance  =  +  1  \distance  +  (T[z  +  l].m  —  T[i  —  l].m) 

9:  end  for 

10:  end  for 

11:  end  for 

12:  end  procedure 


generation  is  then  filled  with  the  best  solution  down  until  the  population  limit  is 
reached. 

The  complexity  of  SPEA2  is  0(MN 2)  due  to  the  sorting  phase  of  the  algorithm. 
This  complexity  more  than  improves  on  the  factorial  complexity  of  the  problem.  The 
only  caveat  to  this  is  how  the  genetic  operators  perform  as,  based  on  their  complexity, 
they  could  dominate  the  algorithm  if  not  constructed  competently.  The  algorithm  is 
described  via  pseudo  code  in  Algorithm  4. 


4-5.2  Strength  Pareto  Evolutionary  Algorithm  II.  The  SPEA2  was  devel¬ 
oped  by  Zitzler  [64]  as  an  improvement  to  the  original  SPEA  algorithm  [65].  SPEA2 
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Figure  4.3:  Crowding  distance  calculation.  Dark  points  are  non-dominated  solutions. 
[14] 
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Algorithm  4  NSGA2 
1:  procedure  NSGA2 
2:  t  \=  0; 

3:  initialize  P( 0)  :=  {ai(0), a^(0)}  E  P1; 

4:  evaluate  P( 0)  :  ($(ai(0)), ...,  $(aM(0))}; 

5:  while  (t(P(t)  7^  true)  do 

6:  P(t)  =  P(t)  U  Q(t) 

7:  T  =  Fast  Nondominated  Sort(i?,i); 

8:  All  non-dominated  fronts  of  Rt 

9:  while  | P'{t)\  <  N  do 

10:  Crowding  Distance  Assignment  {Pi) 

11:  P'(t)  =  P'{t )  U  ^ 

12:  end  while 

13:  Sort(P'(t),  >n) 

14:  P'(t)  —  P^fO  :  N ];  first  N  elements  in  P’(t) 

15:  crossover :  a'k(t)  rSctPc(P(t))] 

16:  mutate :  a'l(t)  := 

17:  evaluate :  P"{t )  :=  (a"(t)), (a"(t))  : 

18:  <F(a,/(t)),...,<h(a"(t)); 

19:  select :  P(t  +  1)  =  s(P"(t )  [J  Q); 

20:  t  :=  t  H-  1; 

21:  end  while 

22:  end  procedure 


uses  a  strength  ranking  procedure  to  stratify  solutions.  Each  solution  is  assigned  a 
strength  value  based  on  the  number  of  solutions  that  dominate  it.  First,  it  is  deter¬ 
mined  how  many  points  dominate  each  solution,  this  is  referred  to  as  the  fitness  value. 
Then  each  dominated  point  is  assigned  the  sum  of  all  the  fitness  values  that  dominate 
it,  this  value  is  then  called  the  strength  value  of  the  solution.  It  is  this  strength  value 
that  is  used  to  rank  the  solution.  The  strength  ranking  procedure  ensures  that  while 
good  solutions  are  kept,  solutions  that  are  more  isolated  (but  still  dominated)  are  also 
kept  in  order  to  ensure  better  exploration  of  the  solution  space. 

After  all  the  solutions  are  ranked,  an  environmental  selection  method  reduces 
the  population  to  a  user  specified  value.  Zitzler  refers  to  this  value  as  the  archive, 
which  is  a  misleading  term.  The  archive  does  not  actually  save  any  information 
from  one  generation  to  another.  Its  purpose  is  to  give  the  user  the  ability  to  control 
how  many  points  should  be  saved  each  generation.  By  contrast,  in  NSGA2  once  the 
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Figure  4.4:  Truncation  procedure  removing  non-dominated  points  [62], 


solutions  are  ranked  the  crowding  distance  method  runs  until  the  next  population  is 
filled.  In  SPEA2,  the  user  specifies  an  archive  size  that  indicates  how  many  points 
should  be  selected  from  the  ranking  procedure  to  be  put  back  into  the  population. 
The  intended  effect  of  this  is  to  reduce  the  number  of,  what  is  hoped  to  be,  redundant 
points  passing  from  one  generation  to  another.  Figure  4.4  illustrates  this  idea  by 
showing  how  non-dominated  points  would  be  removed  from  a  Pareto  front. 

The  run  time  of  the  algorithm  is  dominated  largely  by  the  truncation  operation. 
The  fitness  assignment  procedure  requires  0(N2)  while  the  truncation  operation  is 
0(N3)  where  N  is  the  number  of  individuals.  While  this  is  a  high  level  of  complexity 
it  is  still  preferable  to  the  factorial  complexity  described  in  Section  4.2.  Algorithm  5 
illustrates  the  full  structure  of  SPEA2. 

4-6  VRPTW  Chromosome  Structure 

Any  chromosome  solution  used  in  a  VRP  must  be  able  to  specify  how  many 
vehicles  are  required  and  which  cities  must  be  visited  in  what  order.  The  solution 
chromosome  defines  a  genotype,  which  is  a  code  corresponding  to  a  phenotype  which 
is  the  actual  solution.  In  terms  of  total  information  the  genotype  does  not  need  to 
contain  redundant  or  implied  information.  For  example,  in  the  VRP  it  is  implied 
that  a  route  starts  at  the  depot  and  ends  there.  Encoding  this  information  in  a 
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Algorithm  5  SPEA2 
1:  procedure  SPEA2 

2:  N  :=  populationsize ; 

3:  N  :=  archivesize ; 

4:  t  :  =  0; 

5:  generate  P( 0)  :=  (ai(0), ...,  aM(0)}  G  P‘; 

6:  generate  P( 0)  :=  0; 

7:  evaluateFitness  P(t ); 

8:  evaluateFitness  P(t ); 

9:  while  (t(P(f)  7^  true)  do 

10:  envS election:  P(t  +  1)  =  nonDomPt(P(t),  P(t))-, 

11:  mateS election:  P(t  +  1)  =  Tournament(P(t  +  1)); 

12:  recombine:  P'{t  +  1)  :=  r0r(P(t)); 

13:  mutate:  P"(t  +  1)  :=  mTm(P(t)); 

14:  evaluate:  P(t  +  1)  :  ^{a'[{t)), ...,  <P(a"(t)); 

15:  P(t  +  l)  =  P(t)); 

16:  t  :=  t  F  lj 

17:  end  while 

18:  end  procedure 


chromosome  would  therefore  be  a  waste  of  space.  There  are  three  ways  to  accomplish 
this,  others  could  be  formulated  but  these  have  been  deemed  effective  through  their 
repeated  usage  [39]. 

A  possible  solution  structure  is  a  bit  string  where  every  bit  corresponds  to  an 
edge  2  in  the  solution  (and  every  bit  is  either  one  or  zero  indicating  whether  it  is  or  is 
not  in  the  solution).  This  structure  is  very  simple  but  grows  large  very  quickly  and  the 
organization  requirement  of  the  VRP  lends  itself  more  toward  real  valued  structures 
anyway.  The  second  structure  is  a  single  array  of  real  values  representing  each  target, 
the  order  of  which  indicates  the  order  of  visitation.  Each  route  is  separated  by  zeros 
as  shown  in  Figure  4.5.  This  structure  is  more  efficient  but  still  requires  the  use  of 
separators  to  indicate  where  a  route  begins  and  ends. 

In  [49]  a  structure  for  a  VRP  chromosome  is  defined  that  uses  a  similar  idea  as 
the  array  structure  but  attaches  each  route  to  a  support  structure,  like  that  seen  in 
Figure  4.6.  The  most  beneficial  aspect  of  this  structure  is  that  changes  made  to  a  given 

2An  edge  being  a  connecting  line  between  two  points  in  a  graph. 
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Figure  4.5:  A  possible  chromosome  structure  for  a  VRP. 


route  do  not  require  a  shift  to  the  entire  array  of  values.  In  Tavares  [49],  this  structure 
is  proposed,  and  shown  to  be,  an  affective  structure  especially  for  the  VRPTW.  This 
structure  is  also  used  in  previous  research  by  Slear  [44]  and  Russell  [40]. 

Solution  Genotype  Solution  Phenotype 

Route  1  -*(  3  H  9  M  7  ) 

Route  2  ^  8HioK5) 

Route  3  — ( 4  MAH' 2  "H  6 ) 

Figure  4.6:  GVR  chromosome  structure  for  the  VRP. 


The  GVR  structure  offers  many  attributes  that  make  it  desirable  as  a  chromo¬ 
some  structure.  Its  information  content  does  not  contain  redundancies.  Each  route 
implies  the  existence  of  a  departure  and  return  to  the  depot  even  though  it  is  not 
explicitly  stated.  This  is  made  possible  by  the  support  structure  that  contains  and 
separates  each  route.  It  is  also  desirable  that  infeasible  solutions  are  not  turned  into 
feasible  solutions  by  adding  customers  but  instead  only  by  rearranging  and  removing 
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customers.  The  impact  of  this  is  that  whenever  a  solution  is  checked  for  feasibility 
after  the  addition  of  a  customer  it  can  be  safely  discarded  if  infeasible,  knowing  the 
solution  is  a  dead  end. 


4-1  SRP  Chromosome  Structure 

The  chromosome  structure  used  for  the  SRP  is  essentially  the  same  as  for  the 
VRPTW.  It  consists  of  single  route  definitions  arranged  in  a  support  structure.  The 
only  difference  is  the  arrangement  of  data  within  the  structure.  Since  each  customer 
must  be  visited  by  more  than  one  vehicle  at  a  time  the  SRP  structure  must  also  reflect 
this.  A  diagram  of  this  is  shown  in  Figure  4.7. 

Solution  Genotype  Solution  Phenotype 

Route  1  3  M  5  H  9  ) 

Route  2  XX"  4 X  7 X 9 ) 

Route  3  XX  4  K  5  H  10  ! 

Route  4  ->(  2  H  6  H8  H10! 

Figure  4.7:  Modified  GVR  chromosome  structure  for  the  SRP. 


4-8  VRPTW  Genetic  Operator  Development 

The  operators  described  in  this  section  are  taken  from  several  publications  [33, 
41,50].  The  genetic  operator  alters  a  solution  in  a  random  manner  by  either  randomly 
changing  the  solution  or  optimizing  some  sub-section  of  the  solution.  This  opens  two 
avenues  to  pursue,  simple  operators  applied  many  times  or  the  use  complex  heuristics 
used  to  intelligently  optimize  part  of  a  solution.  Classic  GAs  made  use  of  random 
operators,  however  more  and  more  hybrid  algorithms  incorporate  heuristics  and  local 
search  techniques  to  great  effect  [47].  For  this  research  a  random  crossover  method, 
three  random  mutations,  and  a  heuristic  based  mutation  operator  are  used. 
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4-8.1  Random  Crossover.  Crossover  is  the  genetic  operation  that  occurs 
most  frequently  and  ensures  that  children  created  from  the  process  are  feasible.  The 
operation  takes  two  parents,  a  donor  and  a  receiver.  A  random  selection  of  customers 
is  selected  from  the  donor  and  placed  into  the  first  available  route  in  a  copy  of  the 
receiver,  after  the  customers  in  the  incoming  sub-route  have  been  removed.  The  first 
available  route  is  the  route  that  when  receiving  the  sub-route  does  not  violate  any 
constraints.  If  no  route  exists  then  the  sub-route  is  added  as  a  new  route  after  the 
copy  customers  are  deleted  in  the  receiver.  The  net  result  of  this  is  a  child  that  is  a 
copy  of  the  receiver  but  contains  some  sub  route  section  from  the  donor.  This  process 
is  illustrated  in  Figure  4.8. 


Figure  4.8:  Random  crossover  operator  for  the  VRPTW. 
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4-8.2  Random  Swap  Mutation.  In  swap  mutation  two  random  customers  in 
a  solution  are  swapped  if  doing  so  does  not  violate  constraints.  If  constraints  become 
violated  the  mutation  does  not  proceed.  This  process  is  illustrated  in  Figure  4.9. 


—  Select  Two  Random  Targets 

R1 

12)  (  1  1l)(l3 

V 

R2 

9  10  2  14 

R3 

(5)  15  (  8  )  (  4 

♦ 

R4 

(  6  )(  3 ) ( 16!  (  7  ) 

After  Swap  Mutation 


Figure  4.9:  Random  swap  operator  for  the  VRPTW. 


4-8.3  Random  Inversion  Mutation.  Select  a  random  sub-route  within  a 
solution  and  reverse  the  order  of  visitation.  The  mutation  does  not  proceed  if  this 
results  in  an  invalid  solution.  This  process  is  illustrated  in  Figure  4.10.  The  resultant 
route  may  or  may  not  be  longer  than  the  original  route. 


4-8-4  Random  Insertion  Mutation.  Move  a  random  customer  to  a  random 
location  in  the  solution  while  ensuring  feasibility.  It  is  possible  to  create  a  new  route 
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with  probability  ^  where  V  is  the  number  of  vehicles  [50] .  This  process  is  illustrated 
in  Figure  4.11. 


Select  Random  Target  and 


After  Swap  Mutation 


Figure  4.11:  Random  insertion  operator  for  the  VRPTW. 


4-8.5  Best  Route  Cost  Mutation.  This  operator  randomly  selects  a  route 
within  a  solution  and  optimizes  its  construction  by  rearranging  customers  within  that 
route.  This  is  accomplished  by  first  searching  the  route  and  determining  which  of  the 
customers  is  closest  to  the  depot.  This  customer  then  becomes  the  first  customer.  The 
closest  customer  to  this  customer  that  is  in  the  route  is  then  moved  next  to  the  the 
first  customer  and  so  forth,  creating  a  route  based  on  customer  proximity.  Customers 
that  can  not  be  feasibly  added  are  moved  to  a  new  route,  ft  is  always  assumed  that 
a  single  customer  within  a  route  is  valid,  without  this  assumption  the  problem  would 
not  be  solvable.  The  construction  of  a  new  route  proceeds  in  the  same  way,  placing 
customers  by  order  of  proximity.  This  process  is  illustrated  in  Figure  4.12. 

4-9  SRP  Evolutionary  Operator  Development 

The  SRP  genetic  operators  are  variants  of  the  VRPTW  operators  altered  to 
take  into  account  the  different  structure  of  the  SRP  solutions.  The  difficultly  in 
developing  these  operators  is  ensuring  the  validity  of  the  child  genotype.  Since  the 
chromosome  contains  location  sensitive  information  across  two  dimensions,  as  opposed 
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Randomly  Selected  Route 


After  Mutation 


R3  Length:  127 


R1 


12  1  11  13 


R3 


R4 


R5 


R3  Length:  95 


Figure  4.12:  Best  Route  Cost  mutation  operator. 


to  the  VRPTW  chromosome  which  is  only  sensitive  across  a  single  route,  making  even 
slight  changes  can  cause  invalid  solutions  to  be  created. 


4-9.1  Split  Mutation.  Split  mutation  randomly  selects  a  route  within  the 
SRP  solution  and  attempts  to  reduce  the  total  length  of  that  route  be  eliminating 
unnecessary  target  visitations.  Each  customer  is  satisfied  with  a  certain  number  of 
UAVs  at  its  location,  however  more  can  be  present  than  are  actually  needed.  This  may 
cause  a  route  to  be  longer  than  it  needs  to  be  since  its  divergence  to  an  unnecessary 
target  takes  longer  than  a  direct  route.  The  split  mutation  operation  determines  if  this 
is  occurring  in  a  random  route  and  attempts  to  remove  the  target  from  the  vehicles 
flight  plan.  If  this  operation  then  results  in  an  infeasible  solution  it  is  considered  to 
have  failed,  and  is  not  implemented.  This  process  is  illustrated  in  Figure  4.13. 


4-9.2  Vertical  Swap  Mutation.  The  vertical  swap  operator  swaps  two  differ¬ 
ent  locations  vertically  in  a  given  solution.  This  is  in  contrast  to  the  VRPTW  swap 
mutation  in  section  4.8.2  in  which  the  swapped  targets  can  be  anywhere.  Columns 
within  the  SRP  have  a  close  approximation  to  time  within  the  solution.  It  is  not 
exact  because  distance  information  is  not  contained  within  the  solution,  and  cities  in 
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—  Select  over  satisfied  customer 


After  removal 


Figure  4.13:  Split  mutation  operator  for  the  SRP. 

the  same  column  may  not  actually  be  visited  at  the  same  time.  The  swap  operator 
randomly  selects  a  column  and  two  different  targets  within  that  column.  These  tar¬ 
gets  are  then  swapped  and  feasibility  is  checked.  An  infeasible  solution  is  not  used. 
Figure  4.14  illustrates  this  process.  Note  the  swap  of  customer  four  and  six  between 
vehicle  two  and  four. 

Select  a  column  and 

two  random  locations  After  Mutation 


4-9.3  Random  Crossover  with  Tightening.  The  crossover  operation  must  be 
done  with  particular  care  as  effective  alterations  to  the  solution  are  difficult  to  achieve. 
The  basic  idea  of  the  crossover  operation  is  the  same  as  the  VRPTW  crossover  op¬ 
eration,  from  two  solutions  a  random  route  is  selected  from  each.  This  route  is  then 
added  to  the  other  solution.  The  problem  is,  unlike  the  VRPTW  crossover  operation, 
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subsections  of  a  route  can  not  be  easily  transferred  between  two  solutions.  In  order 
to  compensate  for  this,  the  route  to  be  crossed  is  added  to  the  solution  as  an  entirely 
new  route.  The  solution  then  undergoes  an  operation  called  tightening.  During  this 
operation  the  solution  is  searched  to  determine  what  customers  are  over  satisfied  or 
visited  at  inappropriate  times,  the  customers  in  the  new  route  are  given  preference  for 
staying.  The  resultant  solution  contains  the  additional  information  of  the  crossover 
operation  without  the  redundancy  or  errors  the  operation  would  otherwise  result  in. 
This  process  is  illustrated  in  Figure  4.15. 


4-10  Chapter  Summary 

This  chapter  contains  the  high  level  design  for  the  MOEA  solution  to  the 
VRPTW  and  SRP  models  defined  in  Chapter  III.  An  evolutionary  algorithm  is  chosen 
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due  to  success  in  its  application  from  previous  research  and  analysis  of  the  problem 
of  vehicle  routing.  Two  selection  methods  and  genetic  operators  for  each  problem 
are  devised  and  detailed.  In  the  next  chapter  it  is  shown  how  these  operations  are 
implemented  into  testable  code. 
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V.  Low  Level  Software  Design 


This  chapter  contains  the  low-level  solution  construction  of  a  router  for  the 
VRPTW  and  SRP.  The  design  presented  in  the  previous  chapter  is  mapped 
to  an  implementable  software  design.  The  programming  language,  data  structures, 
and  code  infrastructure  used  are  presented  as  part  of  the  evolution  of  the  research 
conducted. 

5.1  Implementation  Objectives 

The  high  level  GA  design  necessitates  many  elements  of  preplanning  before 
construction  of  a  solution  can  take  place.  These  decisions  include  coding  language 
decisions,  the  use  of  GA  libraries,  and  the  level  of  generalization  required.  The  level 
of  generalization  implies  how  dependent  the  software  is  on  the  problem  and  how 
difficult  it  is  to  transition  the  software  to  different  problems.  The  coded  solutions  for 
the  VRPTW,  SRP,  and  path  planner  consist  of  a  series  of  processing  steps:  read  in 
problem  data,  maintain  customer(target)  information,  construct  individual  solutions, 
apply  the  MOEA  to  those  solutions,  and  return  the  results  in  a  readable  format.  The 
implementation  must  be  able  to  accomplish  these  steps  and  return  results  compatible 
with  the  simulation  software.  The  software  must  adhere  to  object-oriented  standards, 
be  generalized  enough  to  allow  for  the  alteration  of  algorithm  parameters,  and  be 
compatible  with  available  hardware. 

5.2  Selection  of  an  Evolutionary  Computation  Library 

A  basic  structure  exists  within  all  genetic  algorithm  search  techniques.  This 
structure  is  defined  by  the  transition  of  populations  of  solutions  through  a  modification 
and  selection  process.  Due  to  this  structural  concept  it  is  advantageous  to  use  a 
programming  library  or  other  software  utility  that  has  already  been  created  with  this 
basic  structure  in  mind,  hereafter  refereed  to  as  the  infrastructure  of  the  GA.  There 
are  a  variety  of  infrastructure  options  available  for  evolutionary  computation  across 
many  coding  platforms.  Some  of  these  options  are  listed  in  Table  5.1. 
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Table  5.1:  Evolutionary  Algorithm  Infrastructure  Choices. 


Language 

Object 

Oriented 

Advantages 

Disadvantages 

GALib 

C++ 

Yes 

Previous  work  cre¬ 
ated  in  GALib. 
C++  data  structure 
functionality. 

Construction  is  very 
outdated  and  library 
is  no  longer  sup¬ 
ported  by  author. 
Functionality  is  lim¬ 
ited. 

Java  GA 
Package 

Java 

Yes 

Offers  the  portabil¬ 
ity  of  the  Java  en¬ 
gine.  Java  data 

structures  and  pro¬ 
gramming  architec¬ 
ture  easy  work  with. 

Less  efficient  than 
a  C  based  code. 
Not  easy  to  integrate 
with  previous  work 
coded  in  C  code. 

Matlab  GA 
Toolbox 

NA 

No 

Matlab  is  easiest 
language  to  program 
in  and  toolbox  func¬ 
tionality  allows  for 
fast  construction. 

Very  slow  processing 
and  difficult  to  al¬ 
ter  toolbox  function¬ 
ality  past  preset  con¬ 
struction. 

Open 

Beagle 

C++ 

Yes 

C++  data  structure 
are  very  powerful 
and  efficient.  Very 
powerful  functional¬ 
ity  due  to  strict  ob¬ 
ject  oriented  con¬ 
struction.  Library 

is  still  currently  sup¬ 
ported  by  author. 

Object  oriented 

structure  entails 

high  level  of  com¬ 
plexity  requiring 

longer  initial  pro¬ 
duction  time. 

For  this  research  the  Open  Beagle  (OB)  library  was  selected.  Previous  routing 
work  used  the  GALib  library  [41],  however  for  this  research  it  was  determined  that 
a  transition  to  a  more  contemporary  library  would  be  beneficial.  The  OB  library 
contains  very  powerful  and  well  constructed  tools  for  the  creation  of  evolutionary  al¬ 
gorithms.  It  is  written  in  C++  allowing  for  easier  integration  with  existing  simulation 
uses,  all  of  which  are  written  in  C++,  and  the  library  allows  the  use  of  the  vector 
data  structure.  The  library  is  written  in  very  strict  object  oriented  protocol,  meaning 
little  work  is  required  on  the  part  of  the  user  to  get  program  specific  details  inte- 
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grated  into  the  overall  programm  structure,  assuming  they  are  written  to  the  same 
Object  Oriented  (00)  standard.  More  details  concerning  the  implementation  level 
design  aspects  of  OB  can  be  found  in  Appendix  B  and  online  [18].  The  selection  of 
this  infrastructure  drives  the  code  level  requirement  of  all  the  program  components 
as  well  as  the  data  structures  available  (C++  data  structures). 

5. 3  Software  Design 

The  OB  system  is  capable  of  a  variety  of  different  evolutionary  computation 
methods  including  Genetic  Programming  (GP)  and  Evolutionary  Strategies  which 
arises  from  the  object  oriented  structure  of  OB  shown  in  Figure  5.1. 

GA  GP  Other  EC 

Generic  EC  framework 
Object  oriented  foundations 
C++  Standard  Template  Library  (STL) 

X - \ 

Figure  5.1:  Open  Beagle  software  architecture  [17]. 

The  top  level  GA  structure  contains  the  problem  specific  functions  as  seen  in 
Figure  5.2.  Note  the  use  of  different  classes  within  the  OB  structure  to  separate  the 
various  operations.  An  expansion  of  each  of  these  classes  in  the  following  sections 
shows  their  internal  structure. 

5-4  Replacement  Strategy  and  Algorithm  Structure 

Within  OB,  the  algorithm  being  used  is  implemented  as  a  replacement  strategy 
which  selects  individuals  from  the  population  of  children  and  parents  to  be  put  into 
the  next  generation.  As  presented  in  Chapter  IV,  NSGA2  [14]  and  SPEA2  [62]  are 
selected  as  the  replacement  strategies.  Both  algorithms  are  implemented  in  a  similar 
fashion  within  the  overall  OB  structure,  with  the  only  major  difference  being  the 
ranking  and  selection  functions  used.  The  NSGA2  algorithm  is  installed  as  an  option 
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Figure  5.2:  Evolutionary  Algorithm  Structural  Implementation. 


within  the  OB  program  and  was  written  by  Marc  Parizeau.  The  SPEA2  algorithm 
is  based  on  the  coded  algorithm  written  for  the  PISA  system  [5]  and  was  translated 
into  the  OB  form  by  the  author. 

Figure  5.3  describes  both  a  replacement  and  selection  method  contained  within 
the  selection  class.  The  overall  EA  structure  in  Figure  5.2  shows  the  purpose  of  both 
of  these  methods.  The  selection  class  contains  both  methods  within  its  class. 
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Figure  5.3:  Selection  class  structure. 


5-4-1  NSGA2.  The  NSGA2  uses  two  functions  to  rank  and  replace  indi¬ 
viduals:  sortFastND  and  evalCrowdingDistance.  The  sortFastND  function  ranks  all 
solutions  by  non-dominated  rank.  The  sorted  solutions  are  then  put  back  into  in  the 
population  via  the  evalCrowdingDistance  method.  Figure  5.4  shows  the  structural 
implementation  of  the  NSGA2  algorithm. 

The  sorting  function  is  called  on  line  9  and  returns  an  ordered  list  of  individuals. 
The  structure  that  defines  this  list  is  a  vector  of  integer  vectors.  The  first  vector  in 
the  set  contains  the  rank  zero,  non-dominated  members.  The  vectors  that  follow 
contain  the  remaining  dominated  individuals  in  ranking  order  (i.e.  the  second  vector 
contains  rank  one).  Following  the  sort  and  insertion  of  non-dominated  members  into 
the  population  (lines  8  -  15),  the  last  pareto  front  is  inserted  into  the  population 
using  the  crowding  distance  function.  Note  how  changes  to  the  population  are  made 
by  overwriting  existing  members  (line  25  -  27). 

The  algorithm  begins  by  being  passed  the  ioDeme,  which  a  structure  that  con¬ 
tains  populations  of  individuals.  A  copy  is  made  of  the  population  and  all  children 
are  generated  within  the  copy  population  (lines  4-5).  The  sortFastND  function  then 
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void  NSGA20p : : applyAsReplacementStrategy (Deme&  ioDeme,  Context&  ioContext) 

{ 

//  Generate  a  new  generation  of  individuals,  merged  with  the  actual  one. 
Individual: :Bag  lOff springs (ioDeme) ; 
lOff springs  =  generateNextGenerationO ; 

//  Fast  non-dominated  sorting,  followed  by  insertion  of  the  first  Pareto  fronts. 
NSGA20p: : Fronts  IParetoFronts ; 

sortFastND(lParetoFronts ,  ioDeme . size () ,  lOffsprings,  ioContext); 
unsigned  int  HndexDeme=0 ; 

for(unsigned  int  j=0;  j< (IParetoFronts . sizeO -1) ;  ++j)  { 

for(unsigned  int  k=0;  k<lParetoFronts [j] . size () ;  ++k)  { 

ioDeme [lIndexDeme++]  =  lOffsprings  [IParetoFronts  [j]  [k]  ]  ; 

> 

> 


//  Insertion  of  the  last  Pareto  front,  using  crowding  distance 
Individual :: Bag  lLastFrontlndiv; 

for  (unsigned  int  1=0;  KIParetoFronts  .back()  .  sizeO  ;  ++1)  { 

lLastFrontlndiv. push_back(10ff springs [IParetoFronts .back () [1] ] ) ; 

> 

NSGA20p: : Distances  IDistances; 


evalCrowdingDistance(lDistances,  lLastFrontlndiv) ; 
for(unsigned  int  m=0;  HndexDeme< ioDeme . size () ;  ++m)  { 

ioDeme  [HndexDeme++]  =  lLastFrontlndiv  [IDistances  [m]  .  second]  ; 

} 


> 


Figure  5.4:  NSGA2  Code  Structure 

receives  the  new  population  and  a  vector  for  storing  reference  numbers.  The  original 
ioDeme  is  then  modified  by  overwriting  existing  individuals  with  new  ones  from  the 
lOffsprings  structure  based  on  ID  numbers  in  the  IParetoFronts  structure,  updating 
the  population  with  the  first  rank  of  non-dominated  solutions.  The  evalCrowdingDis- 
tance  function  is  then  used  to  evaluate  the  dominated  members  of  the  population. 
The  remaining  population  members  are  updated  with  the  results  from  this  function. 


5-4-2  SPEA2.  The  SPEA2  structure  is  the  same  as  the  NSGA2  structure,  in 
that  the  individuals  are  generated  and  reinserted  in  the  same  manner.  The  difference 
is  how  they  are  selected  for  insertion.  The  code  structure  for  SPEA2  is  shown  in 
Figure  5.5. 

After  a  new  generation  is  created  the  individuals  are  assigned  a  fitness  value  by 
the  calcFitness  function  on  line  11.  Adjacency  matrices  are  then  generated  at  line 
13.  These  matrices  are  used  in  truncation  functions.  The  matrix  data  structure  is 
implemented  as  another  vector  of  vectors  structure.  Based  on  the  number  of  non- 
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void  SPEA20p : : applyAsReplacementStrategy (Deme&  ioDeme,  Context&  ioContext) 

{ 

//  Generate  a  new  generation  of  individuals,  merged  with  the  actual  one. 

Individual: :Bag  lOff springs (ioDeme) ; 
lOff springs  =  generateNextGenerationO ; 

//except  for  the  first  entry,  IDs  of  zero  indicate  the  indiv.  was  removed 
SPEA20p: :vecFitness  UndivIDValues ; 

//calculate  SPEA  fitness  values  for  all  individuals 
this->calcFitness (UndivIDValues ,  lOffsprings,  ioContext); 

//create  the  distance  matrices  used  in  the  truncation  process 
this->calcDistances (lOffsprings . size () ,  lOffsprings) ; 

//truncate  the  population  using  the  strength  and  fitness  environmental  selection  method 
if  (f itness_bucket [0]  >  (*mArchiveSize)  .getWrappedValueO) 

{ 

truncate_nondominated(outFitnessValues ,  inlndividualPool ,  ioContext) ; 

} 

else  if  (inlndividualPool . size ()  >  (*mArchiveSize) .getWrappedValueO) 

{ 

truncate_dominated(outFitnessValues ,  inlndividualPool,  ioContext) ; 

} 

return; 

//add  the  selected  individuals  back  into  the  population 
unsigned  int  row=0; 
unsigned  int  HndexDeme=0 ; 

for(unsigned  int  k=0;  k<lIndivIDValues [0] . size () ;  ++k)  { 

//a  zero  value  indicates  it  was  removed  from  the  population 
if  (UndivIDValues  [0]  [k]  !=  KILL_VALUE) 

{ 

ioDeme  [UndexDeme]  =  lOffsprings  [UndivIDValues  [0]  [k]  ]  ; 

++HndexDeme ; 

}//end  if 

if  (HndexDeme>=ioDeme .  size() )  //stop  if  you  fill  up  the  deme 
break; 

}//end  for 


Figure  5.5:  NSGA2  Code  Structure 

dominated  individuals  and  the  archive  size,  a  process  of  either  truncating  dominated 
or  non-dominated  members  is  pursued  (line  16).  Member  truncation  occurs  via  the 
process  described  in  Section  4.5.2.  The  UndivIDValues  data  structure  initialized  at 
line  8  is  a  vector  which  contains  the  ID  numbers  of  all  members  in  order.  As  individu¬ 
als  are  truncated  the  ID  number  is  then  changed  to  a  negative  one.  When  individuals 
are  inserted  back  into  ioDeme  the  existence  of  this  value  indicates  that  the  value 
should  not  be  included. 
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5.5  VRPTW  Solution  Implementation 

The  VRPTW  solution  uses  an  object  oriented  structure  divided  into  into  3 
classes:  the  customer  information  class,  genetic  operator  class,  and  initialization  class 
as  seen  in  Figure  5.2.  The  selection  methods  are  defined  in  section  5.4. 

The  customer  information  class  contains  information  about  individual  customers 
and  the  methods  for  evaluating  individual  solutions.  The  genetic  operator  class  en¬ 
compasses  all  the  modifying  operations  that  take  place  on  an  individual.  The  ini¬ 
tialization  class  contains  the  method  that  creates  the  first  population.  This  section 
concludes  with  a  description  of  the  heuristic  mutation  procedure;  Best  Cost  Route 
Mutation.  The  remaining  operations  are  described  in  Chapter  IV. 

5.5.1  GVR  Data  Structure.  As  presented  in  Chapter  IV,  Genetic  Vehicle 
Representation  (GVR)  defines  the  structure  of  a  single  solution.  The  requirement 
imposed  by  GVR  is  a  data  structure  that  contains  each  route,  and  within  each  route, 
each  customer,  in  order.  These  requirements  are  met  by  utilizing  the  vector  data 
structure  available  within  the  C++  language,  specifically,  a  vector  of  integer  vectors, 
shown  in  Figure  5.6. 

Data  Structure 

INTEGER  VECTOR 
INTEGER  VECTOR 
INTEGER  VECTOR 
INTEGER  VECTOR 

Figure  5.6:  GVR  Data  Structure  and  Genotype  for  VRPTW. 

The  first  vector  contains  a  set  of  vectors,  each  of  which  represents  a  single 
route.  Each  of  those  vectors  contain  integers  which  represent  identification  numbers 
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for  customers.  The  utility  of  this  method  is  the  use  of  iterators  and  standard  vector 
operations  that  exist  within  C++.  No  customer  information  is  kept  within  an  indi¬ 
vidual,  the  individual  is  only  a  code  of  integers  that  correspond  to  information  kept 
within  the  customerlnfo  class  described  in  the  next  section. 
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Figure  5.7:  Individuals  placement  within  the  OB  system. 


Open  Beagle  stores  each  individual  as  a  sub  class  member  in  a  population  hi¬ 
erarchy  as  shown  in  Figure  5.7.  The  utility  of  this  structure  allows  populations  and 
individuals  to  be  passed  in  a  completely  encased  class  with  the  individual  genotypes 
referenced  as  sub-class  members.  This  not  only  keeps  the  design  in  line  with  software 
design  principles  but  simplifies  the  construction  of  the  GA  by  generalizing  how  the 
solutions  are  passed  between  different  operation. 


5.5.2  Customer  Class.  The  customer  class  contains  solution  evaluation 
methods  and  customer  relevant  data.  The  evaluation  methods  are  called  by  all  levels 
of  the  algorithm  as  every  solution  must  be  checked  for  validity  during  construction 
and  alteration  which  is  why  the  methods  are  contained  in  a  separate  class,  making 
it  accessible  to  all  other  classes.  The  data  structures  of  the  customer  class  contain 
problem  specific  data  used  in  the  evaluation  process  and  for  heuristic  calculations  in 
some  mutation  operator  methods  (i.e.  Best  Cost  Mutation).  Figure  5.8  shows  this 
structural  concept. 
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Figure  5.8:  Customer  class  structure. 


As  problem  data  is  read  in,  customer  objects  are  created.  Each  object  has  prob¬ 
lem  information  associated  with  it  (ID  number,  x  coordinate,  y  coordinate,  demand, 
ready  time,  due  time,  and  service  time).  These  objects  are  stored  within  a  vector  that 
is  kept  by  the  customerlnfo  class.  A  single  customerlnfo  object  is  created  at  startup 
which  is  persistent  throughout  the  run  of  the  GA. 

Two  important  data  structures  are  used  within  the  Customer  class.  These  are 
the  adjacency  matrix  [10]  and  the  proximity  matrix.  Both  are  created  by  defining 
a  vector  of  integer  vectors.  The  adjacency  matrix  contains  distances  between  all 
customer  locations  (i.e.  edge  costs).  The  purpose  of  the  matrix  is  to  serve  as  a 
lookup  table  for  distance  values  during  the  evaluation  procedures,  alleviating  the 
need  to  constantly  recalculate  distances. 

The  proximity  matrix  contains  an  ordered  list  of  the  customers  by  proximity. 
The  row  the  vector  is  in  implies  the  customer  it  is  associated  with.  The  proximity 
matrix  is  used  to  make  the  process  of  rearranging  customers  and  determining  better 
routes  a  more  efficient  process.  Instead  of  repeatedly  searching  through  the  adja¬ 
cency  matrix  to  determine  the  closest  customer  a  single  search  through  a  row  of  the 
proximity  matrix  returns  the  same  information. 


The  customer  class  also  contains  methods  for  the  evaluation  of  solutions  and  the 
verification  of  individual  routes.  These  methods  are  not  only  used  to  determine  the 
path  length  of  an  individual,  but  also  to  validate  that  changes  made  to  a  route  have 
not  resulted  in  an  invalid  solution.  Two  methods  evaluate  entire  solutions;  evaluateS- 
olutioni)  and  validateSolutionQ .  The  evaluation  method  returns  the  actual  values 
associated  with  the  solution  while  the  verification  is  a  boolean  method  that  returns 
true  if  the  solution  is  valid  and  false  if  not.  Each  of  these  methods  employ  route  level 
methods  called  evaluateRoute  and  validateRoute.  The  evaluateRoute  method  deter¬ 
mines  the  total  length  and  wait  time  of  a  route.  validateRoute  is  a  boolean  method 
which  calls  the  evaluateRoute  method  and  ensures  the  values  are  within  problem 
boundaries.  The  solution  evaluation  methods  are  passed  entire  solutions  (vector  of 
integer  vectors)  while  the  route  level  methods  are  passed  only  integer  vectors.  This 
break  up  of  the  methods  is  done  in  accordance  with  software  coding  principles. 

The  design  impact  of  the  customer  class  is  to  keep  all  problem  specific  data 
outside  of  the  actual  genetic  algorithm.  This  reduces  the  amount  of  information  that 
needs  to  be  kept  track  of  by  the  infrastructure.  It  also  ensures  a  generalization  of 
internal  GA  functions  by  requiring  calls  to  the  customerlnfo  class  rather  than  putting 
the  functionality  in  a  single  isolated  GA  function.  New  problem  types  or  evaluation 
procedure  modifications  can  then  be  more  readily  implemented. 

5.5.3  Population  Initialization.  Before  operations  can  be  performed  on  a 
population,  that  population  needs  to  be  created.  The  initialization  process  creates 
a  set  of  random  population  members,  however  the  process  is  conducted  in  a  more 
intelligent  way  than  simply  assigning  random  customers  to  vehicles.  This  process  as 
well  as  its  placement  within  the  initialization  class  is  shown  in  Figure  5.9. 

To  begin,  a  vector  is  assigned  a  random  series  of  numbers  the  size  of  the  number 
of  customers.  Each  number  represents  a  customer  ID.  The  first  individual  in  the  list 
is  added  to  the  first  route.  To  select  the  customer  that  follows,  a  check  is  made  to  the 
proximity  matrix.  The  function  goes  through  this  list  attempting  to  find  the  closest 
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Figure  5.9:  Population  Initialization  class  structure. 


customer  that  can  be  feasibly  added.  When  no  new  customers  can  be  added  via  the 
proximity  matrix  selection  process,  a  new  customer  is  taken  from  the  random  vector 
to  add  to  the  next  route.  It  is  of  course  ensured  that  through  this  process  the  random 
vector  is  updated  as  to  which  members  are  removed.  This  process  is  repeated  until 
all  solutions  for  the  population  are  created. 


5.5.4  Best  Route  Cost  Mutation.  The  best  route  cost  mutation  begins  by 
first  selecting  a  random  route  within  the  solution  to  be  optimized.  Four  vectors  are 
then  created,  a  route  vector,  a  live  copy  vector,  a  proximity  listing  vector,  and  a 
transfer  list  vector.  The  route  vector  stores  the  route  as  it  is  being  constructed.  The 
live  copy  vector  stores  the  previous  route,  its  members  are  deleted  as  they  are  added 
to  the  route  vector.  The  proximity  listing  stores  the  proximity  listing  for  the  customer 
currently  being  searched  for.  The  transfer  list  vector  stores  customers  that  can  not 
be  added  to  the  solution,  these  customers  are  later  added  into  new  routes. 

A  double  nested  for  loop  determines  which  of  the  customers  in  the  route  is  closest 
to  the  depot.  This  becomes  the  first  city  in  the  optimized  route.  An  iterative  process 
then  commences  where  the  search  for  the  next  customer  in  the  route  is  determined  by 
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searching  through  the  proximity  listing  for  that  customer.  A  customer  is  taken  from 
the  proximity  listing  and  added  to  the  route  vector.  If  the  route  is  valid  that  customer 
is  removed  from  the  live  copy  and  the  proximity  listing  is  updated  to  reference  this 
new  customer.  If  a  customer  within  the  live  copy  vector  can  not  be  added  back  into 
the  route  it  is  put  into  the  transfer  list.  If  the  transfer  list  contains  more  than  one 
customer  the  previously  described  process  repeats,  customers  from  the  transfer  vector 
are  added  into  new  routes  based  on  proximity. 


This  process  results  in  an  optimization  of  the  selected  route  and  an  optimal 
placement  of  the  remaining  customers.  The  operator  class  shown  in  Figure  5.10 
contains  only  the  mutation  operator  described  here,  however  all  operators  adhere  to 
this  same  I/O  structure  within  the  operator  class. 


Operator  Class 


L 


Individual  ID 


7 


Figure  5.10:  Operator  class  structure. 


5.6  SRP  Solution  Implementation 

The  SRP  implementation  is  a  modification  of  the  VRPTW  implementation. 
All  of  the  data  structures  within  the  VRPTW  code  are  used  in  the  same  manner  in 
the  SRP  code.  This  decreased  the  production  time  of  the  SRP  solution.  The  SRP 
software  structure  also  remains  the  same,  containing  a  customer  information  class, 
genetic  operator  class,  and  initialization  class.  Each  of  the  methods  within  these 
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classes  are  then  modified  for  use  the  SRP  solution.  The  overall  structure  is  the  same 
as  that  shown  in  Figure  5.2. 

The  customer  information  class  contains  information  about  individual  customers 
and  the  methods  for  evaluating  individual  solutions.  The  genetic  operator  class  en¬ 
compasses  all  the  modifying  operations  that  take  place  on  an  individual.  The  ini¬ 
tialization  class  contains  the  method  that  creates  the  first  population.  This  section 
concludes  with  a  description  of  the  vertical  swap  mutation  procedure  implementation. 

5.6.1  Modified  GVR  Data  Structure.  SRP  uses  the  same  chromosome  struc¬ 
ture  as  the  VRPTW  with  the  difference  that  customers  do  not  appear  only  once.  The 
implementation  remains  a  vector  of  integer  vectors.  Each  integer  vector  represents  a 
single  vehicles  route  within  the  solution  and  each  integer  in  the  route  represents  the 
ID  of  the  customer  to  be  visited. 

Data  Structure 

INTEGER  VECTOR 
INTEGER  VECTOR 
INTEGER  VECTOR 
INTEGER  VECTOR 


Figure  5.11:  Modified  GVR  Data  Structure  and  Genotype  for  SRP. 

5.6.2  Modified  Customer  Class.  The  customer  class  contains  the  customer 
object  and  methods  for  solution  evaluation.  The  information  kept  about  a  single 
customer  is  the  same  as  for  the  VRPTW  with  the  exception  of  customer  demand. 
Since  demand  is  no  longer  satisfied  by  capacity  and  instead  by  the  number  of  UAVs 
on  site,  the  demand  values  meaning  is  changed. 
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Solutions  can  no  longer  be  evaluated  by  route  since  customer  satisfaction  occurs 
across  all  routes  at  once.  Customer  demand  is  also  not  not  satisfied  until  all  vehicles 
arrive,  so  wait  time  is  not  only  based  on  earliest  arrival  but  also  latest  arriving  vehicle 
(wait  time  equals  earliest  arrival  time  minus  latest  arriving  member  time  if  latest 
member  arrives  past  earliest  start  time).  Due  to  this,  evaluation  take  place  by  column 
across  all  routes  at  once. 

The  process  of  evaluation  works  by  first  determining  the  path  length  from  the 
depot  to  the  fist  customer  via  the  adjacency  matrix.  The  cost  value  is  kept  in  a  vector 
the  size  of  the  number  of  vehicles  in  the  solution.  The  method  then  iterates  through 
the  solution  updating  the  cost  vector  for  each  additional  customer.  Wait  times  are 
calculated  by  determining  the  vehicle  arriving  last.  The  method  ends  by  connecting 
all  routes  back  to  the  depot  to  determine  the  final  path  length  and  cost  for  each 
vehicle. 

5. 6.3  Population  Initialization.  The  initial  population  for  the  SRP  is  created 
using  the  same  heuristic  method  as  described  for  the  VRPTW.  A  random  customer 
is  added  to  a  route  followed  by  the  next,  closest  feasible  customer.  The  proximity  is 
determined  via  the  proximity  matrix.  The  modification  for  the  SRP  stems  from  the 
need  to  track  all  routes  as  they  are  being  constructed.  The  process  is  shown  in  Figure 
5.12. 

The  fist  selected  customer  is  added  to  the  number  of  vehicles  needed  to  satisfy 
its  demand.  The  process  then  iteratively  attempts  to  add  new  customers  to  existing 
routes  in  the  solution  until  the  demand  for  that  customer  is  satisfied.  If  the  demand 
is  not  satisfied  after  all  existing  routes  have  been  attempted  a  new  vehicle  is  added 
to  the  solution.  This  process  continues  until  all  customers  have  been  assigned. 

5. 6.4  Vertical  Swap  Mutation.  To  apply  vertical  swap  mutation  two  random 
locations  are  selected  within  the  individual,  ensuring  that  they  occur  at  the  some 
point  in  time  (i.e.  the  same  column).  To  determine  if  a  swap  is  possible  a  copy  of  the 
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Customers  Remaining:  3456789  10  Customers  Remaining  :  5  7  8  9  10 


solution  is  modified  to  include  the  swapped  members.  This  solution  is  then  passed 
to  the  customer  class  for  evaluation.  If  the  solution  returns  as  valid  the  mutation 
is  performed  on  the  actual  solution  and  returned  to  the  population.  If  the  swap  has 
resulted  in  an  invalid  solution  the  mutation  fails,  if  not,  then  the  mutation  is  successful 
and  the  method  terminates. 

5.7  Chapter  Summary 

This  chapter  presented  the  software  solution  design.  This  design  process  entailed 
selection  of  EA  infrastructure  software,  the  separation  of  tasks  within  the  program, 
and  feasible  implementation  designs  for  the  different  genetic  operators.  Some  of  these 
aspects  were  driven  by  the  infrastructure  choice.  The  OB  library  has  a  very  specific 
construction  procedure  for  each  genetic  operator  simplifying  the  implementation  task. 
Outside  the  scope  of  the  infrastructure  is  the  customer  class  which  separates  all  cus- 
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tomer  data  and  evaluation  methods  from  the  main  evolutionary  algorithm.  In  the 
next  chapter  a  series  of  experiments  are  designed  to  validate  these  design  choices. 
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VI.  Experimental  Procedures 


This  chapter  contains  information  about  the  design  of  experiments  performed  as 
part  of  this  thesis.  The  objective  of  the  experiment  design  phase  is  to  develop 
testing  methods  for  the  problem  model  and  solution  design.  The  experiments  divide 
into  three  sections  which  test  the  effectiveness  of  the  routing  software,  the  path  plan¬ 
ner,  and  the  simulator.  The  purpose  of  these  experiments  is  to  validate  the  algorithm 
design  of  the  routing  software  by  applying  standard  VRPTW  benchmark  problems 
and  comparing  the  solutions  to  best  known  solutions.  Modified  versions  of  these  test 
problems  are  also  applied  to  the  SRP  routing  software  in  order  for  the  results  to  be 
comparable. 

6.1  Experimental  Design  Objectives 

The  goal  of  any  experiment  is  to  contribute  evidence  to  the  hypothesis  proposed. 
In  this  case,  there  are  two  hypotheses  to  test  stemming  from  the  objectives  defined 
in  Chapter  I;  the  proposed  solution  design  for  application  to  the  VRPTW  is  valid 
across  a  spectrum  of  benchmark  problems,  and  the  solution  design  for  application  to 
the  SRP  produces  valid  results  comparable  to  those  obtained  in  the  corresponding 
VRPTW  benchmark.  These  objectives  drive  the  experiment  design  such  that  a  set  of 
benchmarks  are  applied  to  the  VRPTW  and  SRP  solutions  resulting  in  a  set  of  valid 
solutions.  These  solutions  then  contain  measurable  metrics  of  total  path  length,  total 
vehicle  count,  total  wait  time,  and  average  path  length  (these  metrics  apply  to  both 
the  VRPTW  and  SRP).  In  the  case  of  the  SRP  the  benchmarks  are  modified  such 
that  customer  demand  is  an  indication  of  vehicle  count  and  not  capacity  demand,  as 
in  the  VRPTW.  Comparison  of  these  metrics  of  performance  allows  for  an  intelligent 
comparison  of  the  solution  process  to  benchmark  problems. 

Accurate  performance  comparisons  require  the  application  of  different  design 
choices  to  the  same  problem,  using  the  same  settings  when  possible.  To  fulfill  this  re¬ 
quirement  genetic  algorithm  settings  are  chosen  and  kept  constant  across  the  spectrum 
of  algorithm  choices.  These  settings  are  determined  through  empirical  experiments 
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deemed  to  best  represent  the  performance  capability  of  the  different  genetic  operators. 
Population  size  and  generation  limit  are  chosen  within  the  desire  to  limit  program 
run  time. 

Three  different  algorithm  designs,  each  with  two  options  for  selection  strategies, 
are  used  in  the  experimental  procedures.  These  three  designs  are  NSGA2,  SPEA2, 
and  a  biased  elitism  algorithm.  The  biased  elitism  algorithm  uses  no  strategy  to 
rank  solutions  instead  using  an  elitist  ordering  procedure  that  is  biased  toward  path 
length.  The  top  number  of  individuals,  equal  to  the  population  size,  are  selected 
from  the  population  after  genetic  alteration.  Each  of  these  designs  is  then  paired 
with  either  a  random  or  tournament  selection  process.  Recall  that  selection  refers  to 
how  solutions  are  selected  for  genetic  mutation.  Tournament  selection  means  some 
number  of  random  individuals  is  selected  from  the  population,  with  replacement, 
and  ranked  (biased  by  path  length)  with  the  top  rank  selected  for  alteration.  The 
SRP  experiments  employ  only  the  use  of  the  tournament  selection  method  as  random 
selection  was  deemed  more  harmful  to  the  SRP  solution  process  from  the  fact  that 
the  genetic  operators  employ  no  local  search  techniques. 

6.2  VRPTW  and  SRP  Experiments 

The  most  commonly  used  benchmarks  for  the  VRPTW  are  the  Solomon  prob¬ 
lems  developed  in  1987  [45].  They  exist  in  three  different  varieties;  a  random  dis¬ 
tribution  of  customers  (R),  clustered  sets  of  customers  (C),  and  hybrid  (RC).  Each 
of  these  three  problems  comes  in  dimensions  of  twenty  five  and  fifty  customers.  In 
order  to  examine  the  effectiveness  of  the  software  as  well  as  the  impact  of  the  multi¬ 
objective  design,  two  problems  from  each  type  are  tested,  listed  in  Table  6.1.  The 
use  of  this  variety  of  problems  illustrates  the  impact  of  problem  type  on  the  solution 
design  as  well  as  solution  performance  in  different  instances.  The  number  designation 
of  each  problem  constitutes  the  time  windows  that  exist  for  that  problem.  Problems 
that  begin  with  a  one,  such  as  R109,  have  small  time  windows,  while  R206  has  much 
larger  time  windows. 
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Table  6.1:  Solomon  test  problem  selections  (Modified  for  SRP). 


Random 

Cluster 

Hybrid 

25  Targets 

R206 

R109 

C103 

C205 

RC107 

RC202 

50  Targets 

R206 

R109 

C103 

C205 

RC107 

RC202 

Each  test  problem  contains  a  set  of  target  coordinates,  target  time  windows, 
and  vehicle  capacity.  The  test  file  is  structured  as  seen  in  Figure  6.1.  The  Euclidean 
distance  between  targets  is  considered  to  be  the  edge  cost.  The  same  problem  selec¬ 
tions  are  applied  to  the  SRP  solution  modified  in  the  demand  column  to  ensure  that 
each  problem  contains  a  realistic  UAV  requirement. 


RC202 

VEHICLE 

NUMBER 

25 

CAPACITY 

1000 

CUSTOMER 

CUST  NO. 

XCOORD. 

YCOORD. 

DEMAND 

READY  TIME 
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0 

40 
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0 

0 
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0 

1 

25 

85 
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2 

22 
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22 
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4 

20 

80 
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20 
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10 

6 

18 
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20 
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7 

15 
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20 

0 
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10 

8 

15 
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10 
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Figure  6.1:  Solomon  test  file  example  of  25  dimension  hybrid  problem. 


Algorithm  effectiveness  varies  greatly  as  different  parameters  within  the  program 
are  tuned.  The  settings  for  each  algorithm  type  were  determined  from  empirical 
analysis  and  literature  review  [33] .  The  operator  percentage  indicates  the  chance  that 
operator  is  used  on  an  individual  during  the  alteration  phase.  The  more  effective 
operators  are  used  more  often  while  the  random  operators  are  used  less.  All  the 
options  for  the  algorithm  used  to  solve  the  VRPTW  problems  are  listed  in  Table  6.2, 
the  option  for  the  SRP  algorithm  are  listed  in  Table  6.3. 
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Table  6.2:  VRPTW  GA  Settings. 


Operator 

Setting 

Random  Crossover 

40% 

Swap  Mutation 

25% 

Inversion  Mutation 

25% 

Insertion  Mutation 

10% 

Best  Route  Cost  Mutation 

40% 

SPEA2  Archive  Size 

80 

Generation  Limit 

1000 

Population  Size 

100 

y  ratio 

2 

The  SRP  software  experiments  use  the  NSGA2  and  biased  elitism  algorithms. 
The  reason  for  this  is  that  results  from  the  VRPTW  reveal  a  consistent  dominance 
of  these  two  methods  over  SPEA2.  Each  algorithm/problem  experiment  is  run  thirty 
times  in  order  to  ensure  reliable  statistical  analysis.  Each  replacement  strategy  uses 
a  tournament  selection  method.  The  population  size  and  operator  application  per¬ 
centages  are  different  from  the  VRPTW  settings  in  order  to  counter  the  SRPs  fragile 
structure.  More  simple  operations  are  performed  to  take  the  place  of  a  few  intelligent 
operations.  Experiments  are  run  against  a  small  subset  of  the  problems  applied  to 
the  VRPTW  (those  entries  bolded  in  Table  6.3). 


Table  6.3:  SRP  GA  Settings. 


Operator 

Setting 

Random  Crossover 

50% 

Split.  Mutation 

25% 

Vertical  Swap  Mutation 

5% 

Generation  Limit 

5000 

Population  Size 

100 

j  ratio 

2 
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6. 3  Testing  Environment 

Experiments  are  performed  on  an  Opteron  248  processor  operating  at  2.2GHz. 
The  system  has  4  GB  off  chip  memory  and  a  128  KB  LI  cache.  The  code  is  compiled 
in  64-bit  using  GNU  C+-I-.  All  test  programs  are  run  using  a  bash  scripting  method 
and  Portable  Batch  System  (PBS)  submission  system. 

6-4  Chapter  Summary 

This  Chapter  justifies  the  experiments  performed  using  the  VRPTW  and  SRP 
design  developed  in  chapter  V.  A  selection  of  VRPTW  problems  representing  the 
spectrum  of  available  benchmarks  are  used  on  both  the  VRPTW  and  SRP,  with  the 
problems  for  the  SRP  modified  slightly  in  terms  of  of  what  the  demand  per  customer 
represents.  I11  the  next  chapter  the  results  from  these  experiments  are  analyzed. 
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VII.  Results  and  Analysis 


This  chapter  contains  results  obtained  from  experimental  procedures  in  Chap¬ 
ter  VI,  analytic  analysis  of  the  results  and  their  contextual  meaning.  The 
VRPTW  section  contains  results  from  experimental  trials  compared  to  literature  re¬ 
sults.  A  statistical  comparison  of  these  results  further  compares  the  performance  of 
the  different  procedures  used.  Results  from  the  SRP  exist  solely  in  the  context  of  this 
investigation  and  are  shown  only  in  terms  of  the  results  obtained.  Statistical  analysis 
further  compares  the  procedures  used. 

7.1  VRPTW  Results 

VRPTW  optimization  occurs  across  three  dimensions  of  total  path  length,  total 
wait  time,  and  number  of  vehicles  used.  Previous  analysis,  and  the  classical  view,  of 
this  problem  attempts  to  optimize  path  length  and  the  number  of  vehicles  used  [33] 
or  path  length  alone  [51].  Experiments  performed  in  this  investigation  do  not  yield 
a  single  solution  optimized  in  any  one  direction  but  rather  a  Pareto  front  of  non- 
dominated  values.  In  order  to  compare  the  results  found  here  to  those  in  the  literature 
they  are  first  shown  in  terms  of  the  best  path  length  found  overall.  The  following  box 
plots  show  the  best  path  lengths  available  at  the  time  of  this  writing  compared  to  a 
distribution  of  values  found  from  experimental  trials  (30  trials).  Each  box  plot  shows 
results  for  a  single  problem  across  six  algorithm  settings:  SPEA2,  NSGA2,  Biased 
Single  Objective;  each  of  which  uses  either  random  or  tournament  selection.  The 
wording  used  to  express  each  of  these  settings  is  shown  in  Table  7.1.  Each  plot  also 
shows  the  best  answer  for  path  length  optimization  found  in  Toth  [51]  and  Diaz  [12]. 

Where  appropriate,  a  Kruskal- Wallis  1  test  is  used  to  further  analyze  perfor¬ 
mance  for  the  different  algorithms  on  specific  problems.  Following  these  box  plots, 
solution  space  plots  are  shown  across  dimensions  of  path  length  and  wait  time  in  order 

1A  Kruskal- Wallis  test  is  a  variance  analysis  tool  used  to  test  the  equality  of  a  population  among 
median  groups.  The  result  of  the  test  is  a  window  over  the  values  of  the  distribution  indicating 
if  other  samples  are  significantly  different.  Here,  Matlab  is  used  to  create  a  visual  of  this  window 
showing  which  sample  are  different  from  each  other  using  an  alpha  of  0.05. 
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Table  7.1:  Box  plot  label  explanations  for  VRPTW  experiments. 


Plot  Definition 

Meaning 

SPEA2  Tourn 

SPEA2  replacement  strategy  using  tournament  selection 
for  genetic  operator  application 

NSGA2  Tourn 

NSGA2  replacement  strategy  using  tournament  selection 
for  genetic  operator  application 

Bias  Single  Tourn 

Biased  Single  Objective  replacement  strategy  using  tour¬ 
nament  selection  for  genetic  operator  application 

SPEA2  Rand 

SPEA2  replacement  strategy  using  random  selection  for 
genetic  operator  application 

NSGA2  Rand 

NSGA2  replacement  strategy  using  random  selection  for 
genetic  operator  application 

Bias  Single  Rand 

Biased  Single  Objective  replacement  strategy  using  random 
selection  for  genetic  operator  application 

to  better  examine  algorithmic  performance.  The  drive  for  this  is  that  observing  only 
path  length  can  be  misleading  when  examining  MOEA  performance.  This  section 
also  divides  into  each  type  of  problem  defined  in  Chapter  VI:  random,  cluster,  and 
hybrid. 

7.1.1  Random  Distribution  Problem.  The  difference  in  performance  between 
high  and  low  dimension  problems  is  readily  observable  by  comparing  Figures  7.1  and 
7.2.  NSGA2  is  observed  to  return  results  closer  to  the  best  answer,  followed  by 
the  biased  single  objective  algorithm,  with  SPEA2  doing  worst.  A  Kruskal-Wallis 
statistical  analysis  performed  in  Matlab  confirms  these  visual  observations  as  seen  in 
Figure  7.3  and  7.4. 

The  performance  observation  per  algorithm  is  repeated  in  the  R206  problem. 
NSGA2  again  manages  to  pull  ahead  in  terms  of  the  path  length  objective  and  along 
with  the  biased  algorithm  approaches  the  best  solution  in  the  higher  dimension  50 
customer  problem  in  Figure  7.6.  Observe  the  consistent  convergence  of  solutions  in 
Figure  7.5  for  NSGA2  using  tournament  selection.  The  performance  of  NSGA2  is 
further  clarified  by  the  statistical  plot  in  Figures  7.7  and  7.8. 
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Figure  7.1:  Trial  results  distribution  for  problem  R109  with  25  customers 


Figure  7.2:  Trial  results  for  random  distribution  problem  R109  with  50  customers 

7.1.2  Cluster  Distribution  Problem.  Within  the  cluster  benchmarks  the 
path  length  objective  becomes  less  consistent  in  returns.  Figure  7.9  shows  all  meth¬ 
ods  closing  in  on  the  best  answer  with  NSGA2  actually  achieving  it  in  a  few  trials. 
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Figure  7.3:  Significance  plot  for  R109  with  Figure  7.4:  Significance  plot  for  R109  with 
25  customers  50  customers 
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Figure  7.5:  Trial  results  for  random  distribution  problem  R206  with  25  customers 


Increasing  the  dimension  of  the  problem,  in  Figure  7.10,  causes  a  return  to  the  per¬ 
formance  seen  so  far,  with  no  algorithm  approaching  the  best  solution.  Statistical 
analysis  shows  only  SPEA2  being  significantly  outperformed  in  Figure  7.11  and  7.12. 

Figure  7.13  shows  convergence  using  the  biased  algorithm  with  a  wide  dispersion 
of  points  using  SPEA2  or  NSGA2.  Figure  7.14  maintains  the  decrease  in  performance 
over  higher  dimensional  problems.  Statistical  analysis  in  Figures  7.15  and  7.16  show 
NSGA2  returning  better  results,  ft  is  interesting  to  note  the  results  of  NSGA2  with 
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Figure  7.6:  Trial  results  for  random  distribution  problem  R206  with  50  customers 

random  selection  returning  with  such  consistent  results.  This  can  be  most  likely  at¬ 
tributed  to  the  nature  of  the  cluster  problems  working  well  with  the  genetic  operators 
used. 


7.1.3  Hybrid  Distribution  Problem.  The  hybrid  problem  would  seem  to 
represent  the  most  difficult  landscape  to  work  in,  however  Figure  7.17  shows  that 
no  algorithm  had  particular  trouble  arriving  at  the  optimal  solution.  This  is  less 
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Figure  7.7:  Significance  plot  for  R206  with  Figure  7.8:  Significance  plot  for  R206  with 
25  customers  50  customers 
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7.9:  Trial  results  for  random  distribution  problem  C103  with  25  customers 
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Figure  7.10:  Trial  results  for  random  distribution  problem  C103  with  50  customers 


true  for  the  fifty  dimension  problem  in  Figure  7.18  as  seen  from  the  results  being 
further  from  the  optimal  value  line.  Statistical  analysis  of  the  results  in  Figure  7.20 
show  NSGA2  and  the  biased  algorithm  returning  consistent  values  with  tournament 
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Figure  7.11:  Significance  plot  for  C103  Figure  7.12:  Significance  plot  for  C103 
with  25  customers  with  50  customers 
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Figure  7.13:  Trial  results  for  random  distribution  problem  C205  with  25  customers 


selection  being  the  deciding  factor  in  superior  performance.  Further  generational 
development  would  most  likely  force  the  solution  closer  to  the  optimal.  Figure  7.19 
shows  NSGA2  achieving  statistically  better  results  by  a  small  margin,  further  leading 
to  the  conclusion  of  its  usefulness  in  developing  solutions.  However,  previous  results 
also  show  consistent  returns  using  the  biased  algorithm,  meaning  no  one  strategy 
dominates  overall. 
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Figure  7.14:  Trial  results  for  random  distribution  problem  C205  with  50  customers 


The  hybrid  problem  in  Figure  7.21  again  shows  convergence  of  the  biased  al¬ 
gorithm  while  NSGA2  and  SPEA2  maintain  a  larger  coverage.  This  is  also  seen  in 
Figure  7.22.  The  comparison  plots  in  Figures  7.23  and  7.24  show  consistent  results 
across  NSGA2  and  the  biased  algorithm.  SPEA2  is  again  beaten  in  this  particular 
performance  measure. 
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Figure  7.15:  Significance  plot  for  C205 
with  25  customers 
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Figure  7.16:  Significance  plot  for  C205 
with  50  customers 
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Figure  7.17:  Trial  results  for  random  distribution  problem  RC107  with  25  customers 


7.2  Impact  of  VRPTW  results 

Results  from  only  observing  the  path  length  objective  can  be  informative  but 
also  slightly  misleading.  It  might  be  assumed  that  SPEA2  is  being  outperformed  in 
all  problem  instances,  and  in  terms  of  path  length  it  is.  However  as  seen  in  Figures 
7.25-7.27  analysis  of  the  non  dominated  front  generated  by  the  NSGA2  and  SPEA2 
trials  shows  a  return  of  results  consistent  with  what  one  would  expect  from  a  multi¬ 
objective  problem. 

The  conclusion  to  be  made  is  that  the  multi-objective  solution  is  effective  in 
returning  a  broad  range  of  results  and  that  these  results  are  pushing  the  front  of  the 
problem.  It  is  therefore  not  odd  that  the  multi-objective  approach  did  not  return  a 
near  optimal  value  for  path  length.  The  returned  value  represents  the  solution  space 
for  the  objectives  selected.  The  fact  that  the  returned  values  are  close  (i.e  within 
10  percent  in  most  cases)  to  the  highest  benchmark  value,  shows  the  validity  of  the 
MOEA  approach  as  being  able  to  find  the  optimal  value  for  a  single  objective,  while 
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Figure  7.18:  Trial  results  for  random  distribution  problem  RC107  with  50  customers 
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Figure  7.19:  Significance  plot  for  RC107 
with  25  customers 
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Figure  7.20:  Significance  plot  for  RC107 
with  50  customers 
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Figure  7.21:  Trial  results  for  random  distribution  problem  RC202  with  25  customers 


Figure  7.22:  Trial  results  for  random  distribution  problem  RC202  with  50  customers 

also  optimizing  across  the  range  of  objectives.  In  short,  it  appears  that  multi- objective 
optimization  is  appropriate  for  this  particular  routing  problem. 
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Figure  7.23:  Significance  plot  for  RC202  Figure  7.24:  Significance  plot  for  RC202 
with  25  customers  with  50  customers 
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Figure  7.25:  Non-dominated  front  comparing  NSGA2  and  SPEA2  for  C205  with  25 
customers 


7. 3  SRP  Results 

SRP  optimization  occurs  across  dimensions  of  total  path  length,  wait  time,  the 
number  of  vehicles  used,  and  average  path  length.  As  this  problem  formulation  is 
unique  to  this  investigation  there  are  no  readily  comparable  results.  Though  the 
problems  differ  in  formulation  the  objectives  remain  the  same  between  the  SRP  and 
VRPTW.  As  such,  the  results  obtained  from  the  VRPTW  solution  are  compared 
in  order  to  illuminate  the  hypothesis  that  the  SRP  represents  a  superior  problem 
model  in  terms  of  individual  vehicle  operation  and  mission  optimization.  As  in  the 
previous  section  results  are  organized  in  box  plots  representing  trial  results  for  the 
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Figure  7.26:  Non-dominated  front  comparing  NSGA2  and  SPEA2  for  R206  with  25 
customers 

selected  problem.  Definitions  for  the  labels  used  in  the  plots  is  found  in  Table  7.2 
A  comparison  of  the  best  results  from  the  VRPTW  solution  are  also  shown  both  in 
terms  of  total  path  length  and  average  path  length  (even  though  average  path  length 
is  not  an  optimized  objective  in  the  VRPTW). 


Table  7.2:  Box  plot  label  explanations  for  SRP  experiments. 


Plot  Definition 

Meaning 

NSGA2  Tourn 

NSGA2  replacement  strategy  using  tournament  selection 
for  genetic  operator  application 

Bias  Single  Tourn 

Biased  Single  Objective  replacement  strategy  using  tour¬ 
nament  selection  for  genetic  operator  application 

7.3.1  Random  Distribution  Problem.  Figure  7.28  compares  the  total  and 
average  path  length  returned  by  NSGA2  and  biased  algorithm.  The  biased  algorithm 
seems  to  be  converging  while  NSGA2  retains  a  larger  spread  of  the  solution  space. 
In  Figure  7.29  the  average  path  lengths  are  compared  to  the  average  path  length 
returns  for  the  VRPTW  using  the  same  algorithm  type.  Neither  NSGA2  or  the  biased 
algorithm  perform  significantly  better  than  one  another,  as  seen  in  Figure  7.30. 
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Figure  7.27:  Non-dominated  front  comparing  NSGA2  and  SPEA2  for  RC205  with  25 
customers 
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Figure  7.28:  Trial  results  for  random  distribution  problem  RC107  with  25  customers 
comparing  total  path  length  and  average  path  length  per  vehicle 


7.3.2  Cluster  Distribution  Problem.  For  the  cluster  distribution  problem, 
Figure  7.31  shows  the  return  of  total  and  average  path  length  between  NSGA2  and 
the  biased  algorithm.  Further  comparison  between  the  same  VRPTW  results  in  Fig¬ 
ure  7.32  shows  a  decreasing  ability  to  handle  this  particular  type  of  problem.  This 
behavior  should  be  expected  as  the  VRPTW  solution  contains  heuristic  operators  that 
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Figure  7.29:  Comparison  of  average  distance  per  vehicle  in  SRP  results  to  VRPTW 
results  for  problem  R109  with  25  customers 
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Figure  7.30:  Statistical  analysis  comparing  SRP  results  to  VRPTW  results  for  prob¬ 
lem  R109  with  25  customers 

deal  specifically  with  clustered  targets,  while  the  SRP  does  not.  Note  in  Figure  ?? 
the  equal  performance  of  the  biased  algorithm  and  NSGA2  for  the  VRPTW.  Even 
with  the  high  constraints  of  the  SRP  problem  model  it  is  still  possible  to  return  per 
vehicle  path  length  of  the  same  distance. 


95 


Total  Path  Length 


Average  Path  Length 


NSGA2Tourn  Bias  Single  Tourn 


Figure  7.31:  Trial  results  for  random  distribution  problem  C107  with  25  customers 
comparing  total  path  length  and  average  path  length  per  vehicle 


Figure  7.32:  Comparison  of  average  distance  per  vehicle  in  SRP  results  to  VRPTW 
results  for  problem  C103  with  25  customers 


7-4  Comparative  and  Extended  Analysis  of  Results 

Results  from  the  VRPTW  experimentation  showed  a  consistent  return  of  results 
across  a  broad  spectrum  of  problems.  Optimization  along  the  path  length  objective 
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Figure  7.33:  Statistical  analysis  comparing  SRP  results  to  VRPTW  results  for  prob¬ 
lem  C103  with  25  customers 

showed  less  than  optimal  results,  however  when  combined  with  a  view  of  the  achieved 
non-dominated  front,  it  is  clear  that  the  MOEA  strategy  is  working  correctly.  Further 
comparison  between  results  shows  the  NSGA2  algorithm  performing  better  or  as  well 
as  the  biased  algorithm.  However,  in  some  cases  (Figure  7.14  and  7.21)  the  biased 
solution  converged  early  while  the  MOEA  approaches  maintained  a  breadth  of  search 
in  the  solution  space.  These  results  lead  to  the  conclusion  that  the  MOEA  solution 
method  developed  and  implemented  here  is  effective  at  optimization  over  a  range  of 
different  problems. 

Even  with  a  multi- objective  design,  optimization  of  the  path  length  objective 
still  approaches  optimal  value.  NSGA2  is  able  to  achieve  a  path  length  value  within 
10  percent  of  the  optimal  value,  and  is  even  closer  in  some  cases  (Figure  7.9).  It  can 
be  concluded  that  even  while  the  algorithm  is  optimizing  across  multiple  objectives 
the  returns  for  a  single  objective  are  no  being  compromised,  as  evidenced  by  NSGA2s 
performance  on  the  various  benchmark  problems. 

With  the  validation  of  the  MOEA  design  in  place,  obtained  through  analysis 
of  results  over  VRPTW  benchmark  problems,  attention  can  then  be  turned  to  to 
the  SRP  problem  model.  Results  for  the  SRP  again  show  consistent  returns  of  total 
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and  average  path  length.  Average  path  length  is  then  compared  to  the  average  path 
length  returns  for  the  VRPTW  solution  in  order  to  show  that  the  SRP  achieves 
comparable  results,  which  it  does.  That  the  average  path  length  returns  for  the  SRP 
are  comparable  indicates  the  merit  of  the  model  as  a  per  vehicle  optimization  strategy. 
The  purpose  of  the  SRP  as  model  is  to  develop  time  constrained  routes  between  many 
different  targets  each  of  which  requires  some  number  of  vehicle  visitations.  It  is  no 
sunrise  that  total  path  length  is  greater  for  the  SRP  returns,  it  would  have  to  be, 
what  is  important  is  that  the  returned  solution  does  not  require  any  one  vehicle  to 
visit  a  large  number  of  points,  as  would  be  the  case  in  the  VRPTW. 

The  benchmarks  used  in  these  experiments  should  also  be  considered  reflective 
of  real  world  problems  and  not  merely  contrived  problems.  A  real  world  mission 
for  a  compliment  of  UAVs  can  conceivably  contain  20  or  more  targets,  to  which 
these  benchmarks  affectively  match.  The  SRP  model  shows  capability  not  only  as 
a  combinatorics  formulation  but  also  as  an  applicable  model  for  real  world  problem 
formulation,  as  the  solutions  shown  here  validate  a  capability  to  return  consistent 
solutions. 

7.5  Chapter  Summary 

This  chapter  presented  the  results  and  analysis  from  the  VRPTW  and  SRP  ex¬ 
periments  defined  in  Chapter  VI.  Results  from  the  VRPTW  experimentation  showed 
consistent  returns  across  the  spectrum  of  problems  applied  to  it  and  results  compa¬ 
rable  to  single  objective  optimization  results  found  in  the  literature.  Drawing  from 
these  results  leads  to  the  conclusion  that  the  VRPTW  algorithmic  solution  design 
constructed  in  this  investigation  represents  a  valid  construct  for  optimizing  VRPTW 
solutions. 

Results  from  the  SRP  experiments  showed  consistent  returns  for  the  two  algo¬ 
rithmic  procedures  examined.  Average  vehicle  length  results  did  not  compare  consis¬ 
tently  favorably  nor  unfavorably  to  the  VRPTW  results.  As  the  objectives  for  both 
problem  models  are  arrived  at  through  different  procedures  it  should  not  be  con- 


eluded  that  the  VRPTW  model  has  outperformed  the  SRP.  Rather,  the  results  show 
the  SRP  model  achieves  per  vehicle  optimization  comparable  to  the  VRPTW  mean¬ 
ing  that  while  the  number  of  vehicles  used  is  greater,  which  would  be  an  expected 
aspect  of  the  models  design,  each  vehicle  does  the  same  amount  of  average  work.  In 
the  next  chapter  the  objectives  for  this  investigation  are  reviewed  and  future  research 
is  discussed. 
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VIII.  Conclusions 


Ruminations  on  the  topics  covered  in  this  research  lead  to  many  different  trains 
of  thought.  This  chapter  summarizes  the  results  obtained  and  analyzed  in  the 
previous  chapter  and  associates  the  objectives  in  the  first  chapter.  A  discussion  of 
future  research  related  to  both  UAV  simulation  and  direct  extensions  of  this  work 
concludes  the  investigation. 

8. 1  Review  of  Accomplishments 

This  research  endeavor  has  yielded  a  number  of  significant  returns: 

•  Development  of  the  Swarm  Routing  Problem  (SRP)  model  is  shown  to  be 
an  effective  and  solvable  problem  model  for  multi  UAV  routing.  The 
SRP  evolves  from  the  initial  definition  of  the  VRPTW  and  is  modified  it  by 
removing  single  target  visitation  restrictions  and  changing  target  sat¬ 
isfaction  to  be  based  on  vehicle  volume  at  a  certain  point  in  time.  These 
changes  cause  the  model  to  treat  the  vehicles  as  members  of  a  sub-swarm 
as  opposed  to  discrete  vehicles.  The  solution  to  this  problem  model  better  rep¬ 
resents  reality  and  offers  an  optimized  allocation  of  UAV  resources. 

•  Implementation  of  MOEA  solution  in  software  using  the  Open  Beagle 
evolutionary  computation  library.  The  benefit  of  this  library  allowed  for 
cross  use  of  many  methods  between  the  VRPTW  and  SRP  solution.  The  soft¬ 
ware  itself  is  also  easily  expandable  and  general  allowing  for  the  testing 
of  a  very  large  number  of  experimentation  settings  and  algorithmic 
domains. 

•  Results  from  the  VRPTW  experiments  show  consistent  returns  often  within 
and  under  10  percent  of  published  results.  The  benefit  of  the  solution 
presented  is  a  consistent  return  of  results  not  constrained  by  problem 
specific  information. 
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•  The  proposed  solution  and  SRP  experimentation  also  shows  the  desired  effect 

of  decreasing  individual  UAV  path  cost  while  still  routing  all  UAVs 

to  all  targets  within  the  time  limits. 

8.1.1  Objective  1:  Develop  SRP  as  new  model  for  UAV  routing.  The  first 
objective,  defined  in  Section  1.4.1,  is  to  modify  the  existing  VRPTW  problem  into 
a  new  type  of  combinatorics  problem  called  the  SRP.  The  drive  from  this  originates 
from  the  VRPTW  being  an  incomplete  model  for  UAV  routing.  This  stems  from 
the  VRPTW  being  a  closer  model  to  truck  routing,  where  each  vehicle  is  seen  as  a 
large  movable  entity  capable  of  visiting  many  targets.  This  is  not  the  case  with  UAV 
routing  which  consists  of  many  small  vehicles  which  satisfy  target  demand  based 
not  on  inherent  capacity  but  volume  on  location  (the  number  of  UAVs  on  site  at  a 
particular  time).  The  VRPTW  is  modified  to  remove  the  restriction  that  each  target 
be  visited  by  only  a  single  vehicle.  The  complexity  of  the  problem  remains  at  least 
as  hard  the  VRPTW. 

8.1.2  Objective  2:  Develop  and  validate  MOEA  solution  to  VRPTW.  The 
second  objective  uses  an  established  combinatorics  problem  to  model  the  routing  of 
UAVs  across  multiple  targets  within  time  constraints.  The  problem  is  discussed  in 
terms  the  multiple  objectives  of  path  length,  vehicle  count,  and  total  wait  time.  A 
MOEA  solution  structure  is  chosen  after  analysis  of  the  problem  concludes  the  very 
fast  growing  solution  space  necessitates  some  form  of  stochastic  solution.  The  MOEA 
solution  is  first  designed  in  terms  of  a  generic  GA  structure  (defining  the  selection  and 
replacement  method)  and  then  a  definition  of  the  specific  operators  used  to  modify 
individual  solutions.  Within  this  design  phase  an  optimization  strategy  involving 
a  local  search  technique  is  created.  This  local  search  takes  place  with  a  mutation 
operator  called  best  route  cost  (see  Section  4.8.5)  where  a  single  route  is  optimized 
across  the  parameter  of  path  length. 

The  use  of  an  MOEA  necessitates  the  use  of  an  evolutionary  algorithm  software 
package,  resulting  in  the  selection  of  the  Open  Beagle  library  written  in  C++.  This 
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library  offers  a  powerful  set  of  tools  for  developing  genetic  algorithms  without  being 
problem  specific.  Genetic  operators  and  an  overall  software  structure  are  designed 
using  this  library  for  the  VRPTW  with  great  success. 

Problems  from  the  Solomon  benchmark  set  are  chosen  in  order  to  validate  the 
solution  and  have  comparable  results  to  other  achievements  in  the  field.  Results  show 
effective  solutions  for  25  and  and  50  dimension  problems  with  large  time  windows  with 
a  degradation  in  performance  for  large  dimension,  highly  constrained  time  window 
problems.  As  benchmark  problems  are  optimized  only  for  the  single  objective  of  path 
length  defining  algorithm  effectiveness  in  terms  of  path  length  is  ineffective,  though 
it  is  used  in  order  to  have  a  comparable  metric.  Resultant  solutions  were  within  10 
percent  of  the  benchmark  value  in  most  cases.  The  net  effect  of  this  experimenta¬ 
tion  shows  the  overall  effectiveness  of  the  MOEA  routing  method  to  develop  “good” 
solutions  to  routing  problems. 

8.1.3  Objective  3:  Develop  and  validate  MOEA  solution  to  SRP.  As  the 
complexity  of  the  SRP  is  at  least  that  of  the  VRPTW,  shown  in  Section  4.2,  a 
stochastic  solution  is  deemed  appropriate.  The  solution  process  parallels  that  of  the 
VRPTW  in  terms  of  algorithm  structure  and  the  code  library  used.  Experiments  used 
the  Solomon  problems  for  the  VRPTW.  These  problems  were  chosen  due  to  there  ease 
of  availability  and  to  have  some  type  of  result  to  compare  to. 

Since  the  structure  of  the  SRP  is  different  from  the  VRPTW  exact  solution 
comparison  is  not  possible.  The  same  objectives  exist  for  both  problems,  however 
since  targets  are  satisfied  in  different  manners  comparing  objectives  of  path  length 
and  vehicle  count  are  not  appropriate.  Instead,  average  path  length  per  vehicle  is 
used  a  comparison  metric.  Results  showed  that  the  proposed  solution  generated  better 
results  over  a  set  number  of  iterations  and  that  average  path  length  was  reduced.  How 
much  better  the  results  are  is  arbitrary  as  the  results  generated  for  these  benchmarks 
is  unique  to  this  work.  The  important  point  to  take  from  this  objective  is  that  by 
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using  the  SRP  problem  model  routing  solutions  can  be  developed  which  “optimize” 
the  path  for  individual  UAVs  while  still  solving  a  large  and  complex  routing  problem. 

8.2  Future  Research  and  Closing  Remarks 

As  mentioned  in  Chapter  II  this  effort  is  one  in  a  long  line  of  research  concerning 
UAV  routing  and  simulation.  With  the  development  of  this  advanced  routing  capa¬ 
bility,  new  developments  in  SO  control  schemes  [32],  and  previous  work  developing 
a  3D  UAV  simulator  the  goal  still  remains  of  developing  a  high  quality  UAV  swarm 
simulator  (operational  in  either  online  or  off  line  mode).  Within  the  confines  of  this 
research  next  step  investigations  within  UAV  routing  would  be: 

1.  Further  research  on  the  SRP  as  a  UAV  routing  problem  model. 

2.  Integration  of  routing  software  with  path  planning  software  [44], 

3.  Full  integration  of  mission  planning  software,  recently  developed  SO  controls 
[32],  and  the  AFIT  simulator  [40].  Each  of  these  pieces  now  exists  and  are 
ready  to  be  combined. 

Outside  the  realm  of  UAV  routing  the  following  topics  could  easily  be  extended 
from  this  work.  These  topics  would  extend  understanding  of  the  VRPTW,  SRP,  and 
MOEA  solution  techniques: 

1.  Development  of  different  heuristic  operators  for  the  VRPTW  in  order  to  further 
optimize  performance.  Many  local  search  techniques  have  already  been  defined 
for  the  VRPTW  such  as  path  savings  strategies  and  problem  relaxation  tech¬ 
niques.  It  would  take  a  small  amount  of  effort  to  convert  one  of  these  search 
techniques  into  a  mutation  operator. 

2.  Alteration  of  design  algorithm  to  handle  other  VRP  variants  to  match  different 
real  world  problem  applications. 
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3.  Expansion  of  MOEA  solution  method  to  handle  very  high  dimension  problems 
(>  100  targets).  Dealing  with  problem  dimensions  higher  than  50  entails  a 
number  of  complexity  and  decision  problems. 

4.  Further  research  on  the  SRP  as  a  combinatorics  problems  to  define  different 
solutions  approaches  and  variants.  Some  possible  solution  approaches  are  de¬ 
scribed  in  Chapter  4.3. 

Work  yet  remains  in  the  field  of  UAV  simulation  and  mission  planning.  The  cur¬ 
rent  pace  of  development  for  experiments  such  as  the  ones  shown  here  be  maintained 
or  increased  if  significant  information  is  to  be  obtained  before  hardware  capable  of 
meeting  the  demands  of  a  real  UAV  swarm  are  developed.  Currently  at  AFIT  research 
within  the  Advanced  Navigation  Technology  (ANT)  Lab  is  focused  on  development 
of  onboard  sensor  and  computer  systems  that  can  provide  a  basis  for  further  swarm 
development.  The  fastest  development  to  a  deploy- able  UAV  swarm  is  achieved  when 
a  full  understanding  and  simulatable  model  of  UAV  swarms  has  been  developed  long 
before  capable  hardware  is. 
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Appendix  A.  Evolutionary  Computation 

Evolutionary  computation  is  a  subfield  within  the  domain  of  natural  computing  [8]. 
Natural  computing  an  area  of  research  in  computer  science  that  attempts  to  emulate 
and  simulate  processes  found  in  nature  for  the  purpose  of  solving  real  world  prob¬ 
lems.  Evolutionary  computation  encompasses  a  wide  number  of  search  techniques. 
Throughout  the  literature  the  term  GA  is  also  used  to  signify  an  evolutionary  algo¬ 
rithm.  In  the  context  of  this  work  both  terms  refer  to  the  same  algorithm  type. 

An  EC  algorithm  uses  a  process  of  manipulation  in  order  to  generate  optimal 
solutions  to  single  or  multi-objective  complex  problems.  This  process  (refereed  to  as 
evolving  a  solution)  involves  using  methods  of  recombination,  mutation,  and  selection 
on  a  population  of  solutions  in  order  to  develop  this  optimal  solution  [2] .  This  process 
is  discussed  in  further  detail  in  section  A.l. 

An  important  aspect  of  this  evolutionary  process  is  the  structure  of  the  chromo¬ 
some  that  defines  an  individual  solution.  This  structure  is  refereed  to  as  a  genotype 
and  is  what  determines  how  the  different  processes  in  the  GA  occur  (the  actual  solu¬ 
tion  is  called  the  phenotype).  The  structure  is  also  what  determines  how  the  solution 
is  represented.  Often  the  structure  of  the  genotype  determines  how  effective  the  GA 
is  at  finding  better  solutions.  In  this  paper  we  present  a  new  type  of  chromosome 
that  possesses  redundant  information  treated  as  a  type  of  memory  for  that  particular 
solution.  By  using  this  structure  the  algorithm  is  able  to  make  more  efficient  use  of 
the  GA  methods.  This  structure  is  incorporated  into  a  GA  and  applied  to  a  Traveling 
Salesman  problem  in  order  to  illustrate  the  effect. 

A.l  Classic  Genetic  Algorithm 

A  standard  GA  is  defined  mathematically  as  an  evolutionary  algorithm  in  [2], 
The  EA  works  via  a  process  of  selection,  crossover,  and  mutation  and  is  based  on  the 
chromosome  structure  the  algorithm  uses.  For  this  paper  we  use  a  TSP  problem  to 
illustrate  these  functions  within  a  GA. 
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A.  1.1  Chromosome  Structure.  The  chromosome  structure  (genotype)  is 
what  determines  the  solution  to  the  problem  being  solved  (phenotype).  The  TSP 
is  a  combinatorics  problem  that  easily  lends  itself  to  being  solved  with  a  Genetic 
Algorithm  due  to  the  way  a  solution  can  be  structured.  There  are  many  variations  to 
the  classic  TSP,  however  in  its  most  general  form  it  is  a  graph  consisting  of  n  vertexes 
and  n-1  edges.  The  goal  is  to  determine  a  path  from  one  vertex,  through  all  the  other 
vertex  points  exactly  once  and  ending  back  at  the  start  point  (called  a  circuit).  The 
objective  of  the  problem  is  to  determine  the  shortest  possible  path  that  accomplishes 
this  goal.  This  is  illustrated  in  Figure  A.l  though  the  graph  edges  are  not  shown.  It 
is  assumed  that  distances  between  cities  is  Euclidean  (i.e.  the  distance  from  1  to  2  is 
less  than  the  distance  from  1  to  5). 


The  chromosome  structure  is  a  code  that  defines  a  particular  solution.  In  [27] 
a  TSP  genotype  is  formatted  as  follows.  All  cities  in  the  graph  are  viewed  as  an 
index.  Each  city  has  a  connection  point  that  indicates  which  city  it  leads  to.  This 
chromosome  structure  as  well  as  the  resultant  TSP  solution  is  illustrated  in  Figure  A. 2. 
This  structure  is  used  to  facilitate  crossover  operations  that  are  discussed  in  the  next 
section. 
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Please  note  that  the  structure  illustrated  is  not  the  exact  structure  that  would 
be  used  in  a  coded  implementation.  Obviously  the  information  indicating  the  index  is 
redundant  and  is  only  shown  for  illustrative  purposes.  The  only  necessary  information 
is  the  second  string  where  the  index  is  implied  by  the  order  of  the  vertexes. 


A.  1.2  Crossover,  Mutation,  and  Selection.  Crossover,  mutation,  and  selec¬ 
tion  are  the  methods  the  GA  uses  that  provide  the  search  functionality.  Each  can  be 
constructed  in  a  variety  of  ways  and  each  impacts  the  effectiveness  of  the  GA  search 
process.  We  continue  to  use  the  TSP  as  the  example  problem. 

Crossover  is  a  process  where  two  solutions  are  ’’crossbred”  resulting  in  two 
children  genotypes,  each  of  which  is  constructed  with  some  aspect  of  its  parent.  This 
process  is  illustrated  in  Figure  A. 3. 

Note  that  the  crossover  operation  is  not  a  static  process.  Often  times  crossover 
results  in  an  invalid  solution.  For  example  when  Parent  II  tries  to  add  its  bottom 
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Figure  A. 2:  A  TSP  Genotype  and  Solution 
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Figure  A. 3:  An  example  of  Crossover 

three  vertexes  to  Child  I  there  is  a  conflict  if  vertex  2  were  to  point  to  vertex  1, 
because  vertex  1  already  points  to  vertex  2.  This  problem  can  be  addressed  in  two 
ways.  Either  the  invalid  solution  can  be  dropped  from  the  population  or  the  solution 
can  be  repaired.  This  repair  result  is  also  shown  in  Figure  A. 3. 

If  only  crossover  operations  occurred  no  new  information  would  ever  be  added 
to  the  system.  Thus  after  a  few  generations  most  of  the  genotypes  would  start  to 
look  the  same  (depending  on  the  selection  process  used).  This  is  the  classic  struggle 
of  exploration  versus  exploitation.  In  order  to  add  new  information  to  the  system  we 
incorporate  the  concept  of  mutation.  Mutation  is  a  random  change  applied  stochasti¬ 
cally  to  a  genotype.  The  result  of  this  operation  is  a  random  tour  that  is  incorporated 
into  the  population  at  a  random  interval.  An  example  of  this  process  is  shown  in  Fig- 
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Figure  A. 4:  An  example  of  Mutation 

ure  A. 4.  Note  that  the  net  effect  of  mutation  is  actually  just  a  switch  of  two  vertex 
points. 

As  new  population  members  are  created  via  Crossover  and  Mutation  a  process 
must  exist  that  eliminates  ineffective  solutions.  In  evolutionary  computation  a  prop¬ 
erty  called  fitness  is  used  to  describe  how  effective  a  solution  is.  In  the  case  of  our 
single  objective  TSP  the  fitness  of  a  solution  is  the  total  length  of  the  circuit.  This 
concept  of  fitness  combined  with  a  deterministic  or  stochastic  selection  method  de¬ 
termines  which  of  the  new  population  members  continue  into  future  generations.  A 
detailed  examination  of  this  method  and  the  others  discussed  in  this  section  can  be 
found  in  [7].  These  operations  are  incorporated  into  a  pseudo  code  description  of  an 
EA  in  Algorithm  6. 
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Algorithm  6  Generic  EA 
1:  procedure  EA 

2:  t  ■=  0; 

3:  generate  P( 0)  :=  (ai(0), a^(0)}  G  JM; 

4:  evaluate  P{ 0)  :  ($(ai(0)), $(aM(0))}; 

5:  while  ( i(P(t )  7^  true )  do 

6:  recombine:  P'{t )  :=  r0r(P(t)); 

7:  mutate :  P"(t )  :=  r®r(P'  (t))\ 

8:  evaluate:  P"(t)  :  $(a"(t)), $(a"(t)); 

9:  select:  P(t  +  1)  =  sQs(P"(t )  [J  Q); 

10:  t  !—  i  H-  1 5 

11:  end  while 

12:  end  procedure 
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Appendix  B.  Implimentation  Documentation 

The  software  library  selected  for  this  investigation  is  the  Open  Beagle  Evolutionary 
Computation  library  [18].  The  library  consists  of  a  complete  evolutionary  compu¬ 
tation  infrastructure  coded  in  object  oriented  C++.  This  library  has  a  number  of 
benefits  over  coding  an  evolutionary  algorithm  for  a  specific  case. 

The  software  is  structured  such  that  each  component  of  an  evolutionary  algo¬ 
rithm  is  represented  within  its  own  class  structure.  Each  structure  is  then  controlled 
and  managed  by  an  overarching  control  structure  that  manages  the  application  of 
each  operation,  the  population,  and  statistical  analysis.  Control  of  all  the  operations 
is  managed  by  a  user  created  XML  file.  This  file  controls  the  structure  of  the  evolu¬ 
tionary  algorithm  to  be  used,  statistical  tools  to  be  applied,  and  how  results  should 
be  stored.  Settings  for  the  algorithm  can  either  be  set  within  the  XML  sheet  or  on 
the  command  line. 

Open  Beagle  also  offers  integrated  printing  of  debug  statements.  Rather  than 
using  standard  output  methods  open  beagle  uses  its  own  structure  for  controlling 
logged  output.  This  allows  for  a  hierarchy  of  debug  information.  By  setting  each 
debug  setting  to  a  level  appropriate  to  its  information  content  total  output  information 
can  then  be  set  from  the  command  line. 

Open  Beagle  offers  an  infinite  expendability  and  modification  capability  within 
the  limits  of  standard  GA  definitions.  By  formulating  each  class  in  only  general  terms 
and  wrapping  the  individuals  within  a  standardized  class  the  implementation  of  any 
GA  design  is  allowed  take  place  without  being  constrained  by  the  Open  Beagle  archi¬ 
tecture.  In  regards  to  the  implementation  discussed  in  this  document  four  important 
aspects  were  created.  The  customer  class,  the  chromosome  data  structure,  and  var¬ 
ious  genetic  operators.  The  construction  concept  of  one  of  the  operators  is  found  in 
Chapter  V. 
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B.l  The  Customer  Class 


The  customer  class  is  created  as  a  persistent  object  at  the  intimidation  phase. 
This  object  then  creates  and  stores  all  the  problem  relevant  data.  Each  operator 
that  requires  the  use  of  a  customer  class  method  is  passed  a  reference  pointer  to  the 
customer  class  object.  In  this  way  all  problem  information  and  evaluation  capability  is 
kept  within  the  customer  class.  This  serves  the  purpose  of  simplifying  the  evaluation 
of  individuals  and  keeping  information  in  a  single  location  as  apposed  to  being  passed 
among  the  population.  Each  population  member  is  only  a  code  of  numbers  whose 
meaning  exists  in  the  customer  class. 

All  methods  within  the  GA  also  call  the  customer  class  in  order  to  evaluate 
solutions.  This  not  only  eliminates  the  need  to  pass  problem  information  between 
classes  it  also  creates  a  single  evaluation  function  that  entire  algorithm  uses  simplifying 
the  construction  process  and  adhering  to  basic  software  engineering  standards. 

B.2  Chromosome  Data  Structure 

The  chromosome  data  structure  had  to  be  created  as  a  standard  one  did  not 
exist  with  the  Open  Beagle  library.  The  design  required  a  data  structure  consisting 
of  a  vector  of  integer  vectors.  To  accomplish  this  the  integer  vector  data  structure 
within  the  Open  Beagle  was  modified  such  that  the  vector  contained  integer  vectors 
instead  of  just  integers.  A  new  print  function  was  also  created  for  displaying  the 
chromosome  in  the  output  hies. 
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